home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / meschach / meschach2.shar < prev    next >
Encoding:
Text File  |  1994-06-07  |  124.0 KB  |  4,848 lines  |  [TEXT/ttxt]
  1. # to unbundle, sh this file (in an empty directory)
  2. echo lufactor.c 1>&2
  3. sed >lufactor.c <<'//GO.SYSIN DD lufactor.c' 's/^-//'
  4. -
  5. -/**************************************************************************
  6. -**
  7. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  8. -**
  9. -**                 Meschach Library
  10. -** 
  11. -** This Meschach Library is provided "as is" without any express 
  12. -** or implied warranty of any kind with respect to this software. 
  13. -** In particular the authors shall not be liable for any direct, 
  14. -** indirect, special, incidental or consequential damages arising 
  15. -** in any way from use of the software.
  16. -** 
  17. -** Everyone is granted permission to copy, modify and redistribute this
  18. -** Meschach Library, provided:
  19. -**  1.  All copies contain this copyright notice.
  20. -**  2.  All modified copies shall carry a notice stating who
  21. -**      made the last modification and the date of such modification.
  22. -**  3.  No charge is made for this software or works derived from it.  
  23. -**      This clause shall not be construed as constraining other software
  24. -**      distributed on the same medium as this software, nor is a
  25. -**      distribution fee considered a charge.
  26. -**
  27. -***************************************************************************/
  28. -
  29. -
  30. -/*
  31. -    Matrix factorisation routines to work with the other matrix files.
  32. -*/
  33. -
  34. -/* LUfactor.c 1.5 11/25/87 */
  35. -static    char    rcsid[] = "$Id: lufactor.c,v 1.7 1994/03/13 23:08:25 des Exp $";
  36. -
  37. -#include    <stdio.h>
  38. -#include    <math.h>
  39. -#include    "matrix.h"
  40. -#include        "matrix2.h"
  41. -
  42. -
  43. -
  44. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  45. -
  46. -/* LUfactor -- gaussian elimination with scaled partial pivoting
  47. -        -- Note: returns LU matrix which is A */
  48. -MAT    *LUfactor(A,pivot)
  49. -MAT    *A;
  50. -PERM    *pivot;
  51. -{
  52. -    u_int    i, j, k, k_max, m, n;
  53. -    int    i_max;
  54. -    Real    **A_v, *A_piv, *A_row;
  55. -    Real  max1, temp;
  56. -    static    VEC    *scale = VNULL;
  57. -
  58. -    if ( A==(MAT *)NULL || pivot==(PERM *)NULL )
  59. -        error(E_NULL,"LUfactor");
  60. -    if ( pivot->size != A->m )
  61. -        error(E_SIZES,"LUfactor");
  62. -    m = A->m;    n = A->n;
  63. -    scale = v_resize(scale,A->m);
  64. -    MEM_STAT_REG(scale,TYPE_VEC);
  65. -    A_v = A->me;
  66. -
  67. -    /* initialise pivot with identity permutation */
  68. -    for ( i=0; i<m; i++ )
  69. -        pivot->pe[i] = i;
  70. -
  71. -    /* set scale parameters */
  72. -    for ( i=0; i<m; i++ )
  73. -    {
  74. -        max1 = 0.0;
  75. -        for ( j=0; j<n; j++ )
  76. -        {
  77. -            temp = fabs(A_v[i][j]);
  78. -            max1 = max(max1,temp);
  79. -        }
  80. -        scale->ve[i] = max1;
  81. -    }
  82. -
  83. -    /* main loop */
  84. -    k_max = min(m,n)-1;
  85. -    for ( k=0; k<k_max; k++ )
  86. -    {
  87. -        /* find best pivot row */
  88. -        max1 = 0.0;    i_max = -1;
  89. -        for ( i=k; i<m; i++ )
  90. -        if ( scale->ve[i] > 0.0 )
  91. -        {
  92. -            temp = fabs(A_v[i][k])/scale->ve[i];
  93. -            if ( temp > max1 )
  94. -            { max1 = temp;    i_max = i;    }
  95. -        }
  96. -        
  97. -        /* if no pivot then ignore column k... */
  98. -        if ( i_max == -1 )
  99. -        continue;
  100. -        
  101. -        /* do we pivot ? */
  102. -        if ( i_max != k )    /* yes we do... */
  103. -        {
  104. -        px_transp(pivot,i_max,k);
  105. -        for ( j=0; j<n; j++ )
  106. -        {
  107. -            temp = A_v[i_max][j];
  108. -            A_v[i_max][j] = A_v[k][j];
  109. -            A_v[k][j] = temp;
  110. -        }
  111. -        }
  112. -        
  113. -        /* row operations */
  114. -        for ( i=k+1; i<m; i++ )    /* for each row do... */
  115. -        {    /* Note: divide by zero should never happen */
  116. -        temp = A_v[i][k] = A_v[i][k]/A_v[k][k];
  117. -        A_piv = &(A_v[k][k+1]);
  118. -        A_row = &(A_v[i][k+1]);
  119. -        if ( k+1 < n )
  120. -            __mltadd__(A_row,A_piv,-temp,(int)(n-(k+1)));
  121. -        /*********************************************
  122. -          for ( j=k+1; j<n; j++ )
  123. -          A_v[i][j] -= temp*A_v[k][j];
  124. -          (*A_row++) -= temp*(*A_piv++);
  125. -          *********************************************/
  126. -        }
  127. -        
  128. -    }
  129. -
  130. -    return A;
  131. -}
  132. -
  133. -
  134. -/* LUsolve -- given an LU factorisation in A, solve Ax=b */
  135. -VEC    *LUsolve(A,pivot,b,x)
  136. -MAT    *A;
  137. -PERM    *pivot;
  138. -VEC    *b,*x;
  139. -{
  140. -    if ( A==(MAT *)NULL || b==(VEC *)NULL || pivot==(PERM *)NULL )
  141. -        error(E_NULL,"LUsolve");
  142. -    if ( A->m != A->n || A->n != b->dim )
  143. -        error(E_SIZES,"LUsolve");
  144. -
  145. -    x = v_resize(x,b->dim);
  146. -    px_vec(pivot,b,x);    /* x := P.b */
  147. -    Lsolve(A,x,x,1.0);    /* implicit diagonal = 1 */
  148. -    Usolve(A,x,x,0.0);    /* explicit diagonal */
  149. -
  150. -    return (x);
  151. -}
  152. -
  153. -/* LUTsolve -- given an LU factorisation in A, solve A^T.x=b */
  154. -VEC    *LUTsolve(LU,pivot,b,x)
  155. -MAT    *LU;
  156. -PERM    *pivot;
  157. -VEC    *b,*x;
  158. -{
  159. -    if ( ! LU || ! b || ! pivot )
  160. -        error(E_NULL,"LUTsolve");
  161. -    if ( LU->m != LU->n || LU->n != b->dim )
  162. -        error(E_SIZES,"LUTsolve");
  163. -
  164. -    x = v_copy(b,x);
  165. -    UTsolve(LU,x,x,0.0);    /* explicit diagonal */
  166. -    LTsolve(LU,x,x,1.0);    /* implicit diagonal = 1 */
  167. -    pxinv_vec(pivot,x,x);    /* x := P^T.tmp */
  168. -
  169. -    return (x);
  170. -}
  171. -
  172. -/* m_inverse -- returns inverse of A, provided A is not too rank deficient
  173. -    -- uses LU factorisation */
  174. -MAT    *m_inverse(A,out)
  175. -MAT    *A, *out;
  176. -{
  177. -    int    i;
  178. -    static VEC    *tmp = VNULL, *tmp2 = VNULL;
  179. -    static MAT    *A_cp = MNULL;
  180. -    static PERM    *pivot = PNULL;
  181. -
  182. -    if ( ! A )
  183. -        error(E_NULL,"m_inverse");
  184. -    if ( A->m != A->n )
  185. -        error(E_SQUARE,"m_inverse");
  186. -    if ( ! out || out->m < A->m || out->n < A->n )
  187. -        out = m_resize(out,A->m,A->n);
  188. -
  189. -    A_cp = m_copy(A,MNULL);
  190. -    tmp = v_resize(tmp,A->m);
  191. -    tmp2 = v_resize(tmp2,A->m);
  192. -    pivot = px_resize(pivot,A->m);
  193. -    MEM_STAT_REG(A_cp,TYPE_MAT);
  194. -    MEM_STAT_REG(tmp, TYPE_VEC);
  195. -    MEM_STAT_REG(tmp2,TYPE_VEC);
  196. -    MEM_STAT_REG(pivot,TYPE_PERM);
  197. -    tracecatch(LUfactor(A_cp,pivot),"m_inverse");
  198. -    for ( i = 0; i < A->n; i++ )
  199. -    {
  200. -        v_zero(tmp);
  201. -        tmp->ve[i] = 1.0;
  202. -        tracecatch(LUsolve(A_cp,pivot,tmp,tmp2),"m_inverse");
  203. -        set_col(out,i,tmp2);
  204. -    }
  205. -
  206. -    return out;
  207. -}
  208. -
  209. -/* LUcondest -- returns an estimate of the condition number of LU given the
  210. -    LU factorisation in compact form */
  211. -double    LUcondest(LU,pivot)
  212. -MAT    *LU;
  213. -PERM    *pivot;
  214. -{
  215. -    static    VEC    *y = VNULL, *z = VNULL;
  216. -    Real    cond_est, L_norm, U_norm, sum;
  217. -    int        i, j, n;
  218. -
  219. -    if ( ! LU || ! pivot )
  220. -    error(E_NULL,"LUcondest");
  221. -    if ( LU->m != LU->n )
  222. -    error(E_SQUARE,"LUcondest");
  223. -    if ( LU->n != pivot->size )
  224. -    error(E_SIZES,"LUcondest");
  225. -
  226. -    n = LU->n;
  227. -    y = v_resize(y,n);
  228. -    z = v_resize(z,n);
  229. -    MEM_STAT_REG(y,TYPE_VEC);
  230. -    MEM_STAT_REG(z,TYPE_VEC);
  231. -
  232. -    for ( i = 0; i < n; i++ )
  233. -    {
  234. -    sum = 0.0;
  235. -    for ( j = 0; j < i; j++ )
  236. -        sum -= LU->me[j][i]*y->ve[j];
  237. -    sum -= (sum < 0.0) ? 1.0 : -1.0;
  238. -    if ( LU->me[i][i] == 0.0 )
  239. -        return HUGE_VAL;
  240. -    y->ve[i] = sum / LU->me[i][i];
  241. -    }
  242. -
  243. -    LTsolve(LU,y,y,1.0);
  244. -    LUsolve(LU,pivot,y,z);
  245. -
  246. -    /* now estimate norm of A (even though it is not directly available) */
  247. -    /* actually computes ||L||_inf.||U||_inf */
  248. -    U_norm = 0.0;
  249. -    for ( i = 0; i < n; i++ )
  250. -    {
  251. -    sum = 0.0;
  252. -    for ( j = i; j < n; j++ )
  253. -        sum += fabs(LU->me[i][j]);
  254. -    if ( sum > U_norm )
  255. -        U_norm = sum;
  256. -    }
  257. -    L_norm = 0.0;
  258. -    for ( i = 0; i < n; i++ )
  259. -    {
  260. -    sum = 1.0;
  261. -    for ( j = 0; j < i; j++ )
  262. -        sum += fabs(LU->me[i][j]);
  263. -    if ( sum > L_norm )
  264. -        L_norm = sum;
  265. -    }
  266. -
  267. -    tracecatch(cond_est = U_norm*L_norm*v_norm_inf(z)/v_norm_inf(y),
  268. -           "LUcondest");
  269. -
  270. -    return cond_est;
  271. -}
  272. //GO.SYSIN DD lufactor.c
  273. echo bkpfacto.c 1>&2
  274. sed >bkpfacto.c <<'//GO.SYSIN DD bkpfacto.c' 's/^-//'
  275. -
  276. -/**************************************************************************
  277. -**
  278. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  279. -**
  280. -**                 Meschach Library
  281. -** 
  282. -** This Meschach Library is provided "as is" without any express 
  283. -** or implied warranty of any kind with respect to this software. 
  284. -** In particular the authors shall not be liable for any direct, 
  285. -** indirect, special, incidental or consequential damages arising 
  286. -** in any way from use of the software.
  287. -** 
  288. -** Everyone is granted permission to copy, modify and redistribute this
  289. -** Meschach Library, provided:
  290. -**  1.  All copies contain this copyright notice.
  291. -**  2.  All modified copies shall carry a notice stating who
  292. -**      made the last modification and the date of such modification.
  293. -**  3.  No charge is made for this software or works derived from it.  
  294. -**      This clause shall not be construed as constraining other software
  295. -**      distributed on the same medium as this software, nor is a
  296. -**      distribution fee considered a charge.
  297. -**
  298. -***************************************************************************/
  299. -
  300. -
  301. -/*
  302. -    Matrix factorisation routines to work with the other matrix files.
  303. -*/
  304. -
  305. -static    char    rcsid[] = "$Id: bkpfacto.c,v 1.7 1994/01/13 05:45:50 des Exp $";
  306. -
  307. -#include    <stdio.h>
  308. -#include    <math.h>
  309. -#include    "matrix.h"
  310. -#include        "matrix2.h"
  311. -
  312. -#define    btos(x)    ((x) ? "TRUE" : "FALSE")
  313. -
  314. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  315. -
  316. -#define alpha    0.6403882032022076 /* = (1+sqrt(17))/8 */
  317. -
  318. -/* sqr -- returns square of x -- utility function */
  319. -double    sqr(x)
  320. -double    x;
  321. -{    return x*x;    }
  322. -
  323. -/* interchange -- a row/column swap routine */
  324. -static void interchange(A,i,j)
  325. -MAT    *A;    /* assumed != NULL & also SQUARE */
  326. -int    i, j;    /* assumed in range */
  327. -{
  328. -    Real    **A_me, tmp;
  329. -    int    k, n;
  330. -
  331. -    A_me = A->me;    n = A->n;
  332. -    if ( i == j )
  333. -        return;
  334. -    if ( i > j )
  335. -    {    k = i;    i = j;    j = k;    }
  336. -    for ( k = 0; k < i; k++ )
  337. -    {
  338. -        /* tmp = A_me[k][i]; */
  339. -        tmp = m_entry(A,k,i);
  340. -        /* A_me[k][i] = A_me[k][j]; */
  341. -        m_set_val(A,k,i,m_entry(A,k,j));
  342. -        /* A_me[k][j] = tmp; */
  343. -        m_set_val(A,k,j,tmp);
  344. -    }
  345. -    for ( k = j+1; k < n; k++ )
  346. -    {
  347. -        /* tmp = A_me[j][k]; */
  348. -        tmp = m_entry(A,j,k);
  349. -        /* A_me[j][k] = A_me[i][k]; */
  350. -        m_set_val(A,j,k,m_entry(A,i,k));
  351. -        /* A_me[i][k] = tmp; */
  352. -        m_set_val(A,i,k,tmp);
  353. -    }
  354. -    for ( k = i+1; k < j; k++ )
  355. -    {
  356. -        /* tmp = A_me[k][j]; */
  357. -        tmp = m_entry(A,k,j);
  358. -        /* A_me[k][j] = A_me[i][k]; */
  359. -        m_set_val(A,k,j,m_entry(A,i,k));
  360. -        /* A_me[i][k] = tmp; */
  361. -        m_set_val(A,i,k,tmp);
  362. -    }
  363. -    /* tmp = A_me[i][i]; */
  364. -    tmp = m_entry(A,i,i);
  365. -    /* A_me[i][i] = A_me[j][j]; */
  366. -    m_set_val(A,i,i,m_entry(A,j,j));
  367. -    /* A_me[j][j] = tmp; */
  368. -    m_set_val(A,j,j,tmp);
  369. -}
  370. -
  371. -/* BKPfactor -- Bunch-Kaufman-Parlett factorisation of A in-situ
  372. -    -- A is factored into the form P'AP = MDM' where 
  373. -    P is a permutation matrix, M lower triangular and D is block
  374. -    diagonal with blocks of size 1 or 2
  375. -    -- P is stored in pivot; blocks[i]==i iff D[i][i] is a block */
  376. -MAT    *BKPfactor(A,pivot,blocks)
  377. -MAT    *A;
  378. -PERM    *pivot, *blocks;
  379. -{
  380. -    int    i, j, k, n, onebyone, r;
  381. -    Real    **A_me, aii, aip1, aip1i, lambda, sigma, tmp;
  382. -    Real    det, s, t;
  383. -
  384. -    if ( ! A || ! pivot || ! blocks )
  385. -        error(E_NULL,"BKPfactor");
  386. -    if ( A->m != A->n )
  387. -        error(E_SQUARE,"BKPfactor");
  388. -    if ( A->m != pivot->size || pivot->size != blocks->size )
  389. -        error(E_SIZES,"BKPfactor");
  390. -
  391. -    n = A->n;
  392. -    A_me = A->me;
  393. -    px_ident(pivot);    px_ident(blocks);
  394. -
  395. -    for ( i = 0; i < n; i = onebyone ? i+1 : i+2 )
  396. -    {
  397. -        /* printf("# Stage: %d\n",i); */
  398. -        aii = fabs(m_entry(A,i,i));
  399. -        lambda = 0.0;    r = (i+1 < n) ? i+1 : i;
  400. -        for ( k = i+1; k < n; k++ )
  401. -        {
  402. -            tmp = fabs(m_entry(A,i,k));
  403. -            if ( tmp >= lambda )
  404. -            {
  405. -            lambda = tmp;
  406. -            r = k;
  407. -            }
  408. -        }
  409. -        /* printf("# lambda = %g, r = %d\n", lambda, r); */
  410. -        /* printf("# |A[%d][%d]| = %g\n",r,r,fabs(m_entry(A,r,r))); */
  411. -
  412. -        /* determine if 1x1 or 2x2 block, and do pivoting if needed */
  413. -        if ( aii >= alpha*lambda )
  414. -        {
  415. -            onebyone = TRUE;
  416. -            goto dopivot;
  417. -        }
  418. -        /* compute sigma */
  419. -        sigma = 0.0;
  420. -        for ( k = i; k < n; k++ )
  421. -        {
  422. -            if ( k == r )
  423. -            continue;
  424. -            tmp = ( k > r ) ? fabs(m_entry(A,r,k)) :
  425. -            fabs(m_entry(A,k,r));
  426. -            if ( tmp > sigma )
  427. -            sigma = tmp;
  428. -        }
  429. -        if ( aii*sigma >= alpha*sqr(lambda) )
  430. -            onebyone = TRUE;
  431. -        else if ( fabs(m_entry(A,r,r)) >= alpha*sigma )
  432. -        {
  433. -            /* printf("# Swapping rows/cols %d and %d\n",i,r); */
  434. -            interchange(A,i,r);
  435. -            px_transp(pivot,i,r);
  436. -            onebyone = TRUE;
  437. -        }
  438. -        else
  439. -        {
  440. -            /* printf("# Swapping rows/cols %d and %d\n",i+1,r); */
  441. -            interchange(A,i+1,r);
  442. -            px_transp(pivot,i+1,r);
  443. -            px_transp(blocks,i,i+1);
  444. -            onebyone = FALSE;
  445. -        }
  446. -        /* printf("onebyone = %s\n",btos(onebyone)); */
  447. -        /* printf("# Matrix so far (@checkpoint A) =\n"); */
  448. -        /* m_output(A); */
  449. -        /* printf("# pivot =\n");    px_output(pivot); */
  450. -        /* printf("# blocks =\n");    px_output(blocks); */
  451. -
  452. -dopivot:
  453. -        if ( onebyone )
  454. -        {   /* do one by one block */
  455. -            if ( m_entry(A,i,i) != 0.0 )
  456. -            {
  457. -            aii = m_entry(A,i,i);
  458. -            for ( j = i+1; j < n; j++ )
  459. -            {
  460. -                tmp = m_entry(A,i,j)/aii;
  461. -                for ( k = j; k < n; k++ )
  462. -                m_sub_val(A,j,k,tmp*m_entry(A,i,k));
  463. -                m_set_val(A,i,j,tmp);
  464. -            }
  465. -            }
  466. -        }
  467. -        else /* onebyone == FALSE */
  468. -        {   /* do two by two block */
  469. -            det = m_entry(A,i,i)*m_entry(A,i+1,i+1)-sqr(m_entry(A,i,i+1));
  470. -            /* Must have det < 0 */
  471. -            /* printf("# det = %g\n",det); */
  472. -            aip1i = m_entry(A,i,i+1)/det;
  473. -            aii = m_entry(A,i,i)/det;
  474. -            aip1 = m_entry(A,i+1,i+1)/det;
  475. -            for ( j = i+2; j < n; j++ )
  476. -            {
  477. -            s = - aip1i*m_entry(A,i+1,j) + aip1*m_entry(A,i,j);
  478. -            t = - aip1i*m_entry(A,i,j) + aii*m_entry(A,i+1,j);
  479. -            for ( k = j; k < n; k++ )
  480. -                m_sub_val(A,j,k,m_entry(A,i,k)*s + m_entry(A,i+1,k)*t);
  481. -            m_set_val(A,i,j,s);
  482. -            m_set_val(A,i+1,j,t);
  483. -            }
  484. -        }
  485. -        /* printf("# Matrix so far (@checkpoint B) =\n"); */
  486. -        /* m_output(A); */
  487. -        /* printf("# pivot =\n");    px_output(pivot); */
  488. -        /* printf("# blocks =\n");    px_output(blocks); */
  489. -    }
  490. -
  491. -    /* set lower triangular half */
  492. -    for ( i = 0; i < A->m; i++ )
  493. -        for ( j = 0; j < i; j++ )
  494. -        m_set_val(A,i,j,m_entry(A,j,i));
  495. -
  496. -    return A;
  497. -}
  498. -
  499. -/* BKPsolve -- solves A.x = b where A has been factored a la BKPfactor()
  500. -    -- returns x, which is created if NULL */
  501. -VEC    *BKPsolve(A,pivot,block,b,x)
  502. -MAT    *A;
  503. -PERM    *pivot, *block;
  504. -VEC    *b, *x;
  505. -{
  506. -    static VEC    *tmp=VNULL;    /* dummy storage needed */
  507. -    int    i, j, n, onebyone;
  508. -    Real    **A_me, a11, a12, a22, b1, b2, det, sum, *tmp_ve, tmp_diag;
  509. -
  510. -    if ( ! A || ! pivot || ! block || ! b )
  511. -        error(E_NULL,"BKPsolve");
  512. -    if ( A->m != A->n )
  513. -        error(E_SQUARE,"BKPsolve");
  514. -    n = A->n;
  515. -    if ( b->dim != n || pivot->size != n || block->size != n )
  516. -        error(E_SIZES,"BKPsolve");
  517. -    x = v_resize(x,n);
  518. -    tmp = v_resize(tmp,n);
  519. -    MEM_STAT_REG(tmp,TYPE_VEC);
  520. -
  521. -    A_me = A->me;    tmp_ve = tmp->ve;
  522. -
  523. -    px_vec(pivot,b,tmp);
  524. -    /* solve for lower triangular part */
  525. -    for ( i = 0; i < n; i++ )
  526. -    {
  527. -        sum = v_entry(tmp,i);
  528. -        if ( block->pe[i] < i )
  529. -            for ( j = 0; j < i-1; j++ )
  530. -            sum -= m_entry(A,i,j)*v_entry(tmp,j);
  531. -        else
  532. -            for ( j = 0; j < i; j++ )
  533. -            sum -= m_entry(A,i,j)*v_entry(tmp,j);
  534. -        v_set_val(tmp,i,sum);
  535. -    }
  536. -    /* printf("# BKPsolve: solving L part: tmp =\n");    v_output(tmp); */
  537. -    /* solve for diagonal part */
  538. -    for ( i = 0; i < n; i = onebyone ? i+1 : i+2 )
  539. -    {
  540. -        onebyone = ( block->pe[i] == i );
  541. -        if ( onebyone )
  542. -        {
  543. -            tmp_diag = m_entry(A,i,i);
  544. -            if ( tmp_diag == 0.0 )
  545. -            error(E_SING,"BKPsolve");
  546. -            /* tmp_ve[i] /= tmp_diag; */
  547. -            v_set_val(tmp,i,v_entry(tmp,i) / tmp_diag);
  548. -        }
  549. -        else
  550. -        {
  551. -            a11 = m_entry(A,i,i);
  552. -            a22 = m_entry(A,i+1,i+1);
  553. -            a12 = m_entry(A,i+1,i);
  554. -            b1 = v_entry(tmp,i);    b2 = v_entry(tmp,i+1);
  555. -            det = a11*a22-a12*a12;    /* < 0 : see BKPfactor() */
  556. -            if ( det == 0.0 )
  557. -            error(E_SING,"BKPsolve");
  558. -            det = 1/det;
  559. -            v_set_val(tmp,i,det*(a22*b1-a12*b2));
  560. -            v_set_val(tmp,i+1,det*(a11*b2-a12*b1));
  561. -        }
  562. -    }
  563. -    /* printf("# BKPsolve: solving D part: tmp =\n");    v_output(tmp); */
  564. -    /* solve for transpose of lower traingular part */
  565. -    for ( i = n-1; i >= 0; i-- )
  566. -    {    /* use symmetry of factored form to get stride 1 */
  567. -        sum = v_entry(tmp,i);
  568. -        if ( block->pe[i] > i )
  569. -            for ( j = i+2; j < n; j++ )
  570. -            sum -= m_entry(A,i,j)*v_entry(tmp,j);
  571. -        else
  572. -            for ( j = i+1; j < n; j++ )
  573. -            sum -= m_entry(A,i,j)*v_entry(tmp,j);
  574. -        v_set_val(tmp,i,sum);
  575. -    }
  576. -
  577. -    /* printf("# BKPsolve: solving L^T part: tmp =\n");v_output(tmp); */
  578. -    /* and do final permutation */
  579. -    x = pxinv_vec(pivot,tmp,x);
  580. -
  581. -    return x;
  582. -}
  583. -
  584. -        
  585. -
  586. //GO.SYSIN DD bkpfacto.c
  587. echo chfactor.c 1>&2
  588. sed >chfactor.c <<'//GO.SYSIN DD chfactor.c' 's/^-//'
  589. -
  590. -/**************************************************************************
  591. -**
  592. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  593. -**
  594. -**                 Meschach Library
  595. -** 
  596. -** This Meschach Library is provided "as is" without any express 
  597. -** or implied warranty of any kind with respect to this software. 
  598. -** In particular the authors shall not be liable for any direct, 
  599. -** indirect, special, incidental or consequential damages arising 
  600. -** in any way from use of the software.
  601. -** 
  602. -** Everyone is granted permission to copy, modify and redistribute this
  603. -** Meschach Library, provided:
  604. -**  1.  All copies contain this copyright notice.
  605. -**  2.  All modified copies shall carry a notice stating who
  606. -**      made the last modification and the date of such modification.
  607. -**  3.  No charge is made for this software or works derived from it.  
  608. -**      This clause shall not be construed as constraining other software
  609. -**      distributed on the same medium as this software, nor is a
  610. -**      distribution fee considered a charge.
  611. -**
  612. -***************************************************************************/
  613. -
  614. -
  615. -/*
  616. -    Matrix factorisation routines to work with the other matrix files.
  617. -*/
  618. -
  619. -/* CHfactor.c 1.2 11/25/87 */
  620. -static    char    rcsid[] = "$Id: chfactor.c,v 1.2 1994/01/13 05:36:36 des Exp $";
  621. -
  622. -#include    <stdio.h>
  623. -#include    <math.h>
  624. -#include    "matrix.h"
  625. -#include        "matrix2.h"
  626. -
  627. -
  628. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  629. -
  630. -/* CHfactor -- Cholesky L.L' factorisation of A in-situ */
  631. -MAT    *CHfactor(A)
  632. -MAT    *A;
  633. -{
  634. -    u_int    i, j, k, n;
  635. -    Real    **A_ent, *A_piv, *A_row, sum, tmp;
  636. -
  637. -    if ( A==(MAT *)NULL )
  638. -        error(E_NULL,"CHfactor");
  639. -    if ( A->m != A->n )
  640. -        error(E_SQUARE,"CHfactor");
  641. -    n = A->n;    A_ent = A->me;
  642. -
  643. -    for ( k=0; k<n; k++ )
  644. -    {    
  645. -        /* do diagonal element */
  646. -        sum = A_ent[k][k];
  647. -        A_piv = A_ent[k];
  648. -        for ( j=0; j<k; j++ )
  649. -        {
  650. -            /* tmp = A_ent[k][j]; */
  651. -            tmp = *A_piv++;
  652. -            sum -= tmp*tmp;
  653. -        }
  654. -        if ( sum <= 0.0 )
  655. -            error(E_POSDEF,"CHfactor");
  656. -        A_ent[k][k] = sqrt(sum);
  657. -
  658. -        /* set values of column k */
  659. -        for ( i=k+1; i<n; i++ )
  660. -        {
  661. -            sum = A_ent[i][k];
  662. -            A_piv = A_ent[k];
  663. -            A_row = A_ent[i];
  664. -            sum -= __ip__(A_row,A_piv,(int)k);
  665. -            /************************************************
  666. -            for ( j=0; j<k; j++ )
  667. -                sum -= A_ent[i][j]*A_ent[k][j];
  668. -                sum -= (*A_row++)*(*A_piv++);
  669. -            ************************************************/
  670. -            A_ent[j][i] = A_ent[i][j] = sum/A_ent[k][k];
  671. -        }
  672. -    }
  673. -
  674. -    return (A);
  675. -}
  676. -
  677. -
  678. -/* CHsolve -- given a CHolesky factorisation in A, solve A.x=b */
  679. -VEC    *CHsolve(A,b,x)
  680. -MAT    *A;
  681. -VEC    *b,*x;
  682. -{
  683. -    if ( A==(MAT *)NULL || b==(VEC *)NULL )
  684. -        error(E_NULL,"CHsolve");
  685. -    if ( A->m != A->n || A->n != b->dim )
  686. -        error(E_SIZES,"CHsolve");
  687. -    x = v_resize(x,b->dim);
  688. -    Lsolve(A,b,x,0.0);
  689. -    Usolve(A,x,x,0.0);
  690. -
  691. -    return (x);
  692. -}
  693. -
  694. -/* LDLfactor -- L.D.L' factorisation of A in-situ */
  695. -MAT    *LDLfactor(A)
  696. -MAT    *A;
  697. -{
  698. -    u_int    i, k, n, p;
  699. -    Real    **A_ent;
  700. -    Real d, sum;
  701. -    static VEC    *r = VNULL;
  702. -
  703. -    if ( ! A )
  704. -        error(E_NULL,"LDLfactor");
  705. -    if ( A->m != A->n )
  706. -        error(E_SQUARE,"LDLfactor");
  707. -    n = A->n;    A_ent = A->me;
  708. -    r = v_resize(r,n);
  709. -    MEM_STAT_REG(r,TYPE_VEC);
  710. -
  711. -    for ( k = 0; k < n; k++ )
  712. -    {
  713. -        sum = 0.0;
  714. -        for ( p = 0; p < k; p++ )
  715. -        {
  716. -            r->ve[p] = A_ent[p][p]*A_ent[k][p];
  717. -            sum += r->ve[p]*A_ent[k][p];
  718. -        }
  719. -        d = A_ent[k][k] -= sum;
  720. -
  721. -        if ( d == 0.0 )
  722. -            error(E_SING,"LDLfactor");
  723. -        for ( i = k+1; i < n; i++ )
  724. -        {
  725. -            sum = __ip__(A_ent[i],r->ve,(int)k);
  726. -            /****************************************
  727. -            sum = 0.0;
  728. -            for ( p = 0; p < k; p++ )
  729. -            sum += A_ent[i][p]*r->ve[p];
  730. -            ****************************************/
  731. -            A_ent[i][k] = (A_ent[i][k] - sum)/d;
  732. -        }
  733. -    }
  734. -
  735. -    return A;
  736. -}
  737. -
  738. -VEC    *LDLsolve(LDL,b,x)
  739. -MAT    *LDL;
  740. -VEC    *b, *x;
  741. -{
  742. -    if ( ! LDL || ! b )
  743. -        error(E_NULL,"LDLsolve");
  744. -    if ( LDL->m != LDL->n )
  745. -        error(E_SQUARE,"LDLsolve");
  746. -    if ( LDL->m != b->dim )
  747. -        error(E_SIZES,"LDLsolve");
  748. -    x = v_resize(x,b->dim);
  749. -
  750. -    Lsolve(LDL,b,x,1.0);
  751. -    Dsolve(LDL,x,x);
  752. -    LTsolve(LDL,x,x,1.0);
  753. -
  754. -    return x;
  755. -}
  756. -
  757. -/* MCHfactor -- Modified Cholesky L.L' factorisation of A in-situ */
  758. -MAT    *MCHfactor(A,tol)
  759. -MAT    *A;
  760. -double  tol;
  761. -{
  762. -    u_int    i, j, k, n;
  763. -    Real    **A_ent, *A_piv, *A_row, sum, tmp;
  764. -
  765. -    if ( A==(MAT *)NULL )
  766. -        error(E_NULL,"MCHfactor");
  767. -    if ( A->m != A->n )
  768. -        error(E_SQUARE,"MCHfactor");
  769. -    if ( tol <= 0.0 )
  770. -            error(E_RANGE,"MCHfactor");
  771. -    n = A->n;    A_ent = A->me;
  772. -
  773. -    for ( k=0; k<n; k++ )
  774. -    {    
  775. -        /* do diagonal element */
  776. -        sum = A_ent[k][k];
  777. -        A_piv = A_ent[k];
  778. -        for ( j=0; j<k; j++ )
  779. -        {
  780. -            /* tmp = A_ent[k][j]; */
  781. -            tmp = *A_piv++;
  782. -            sum -= tmp*tmp;
  783. -        }
  784. -        if ( sum <= tol )
  785. -            sum = tol;
  786. -        A_ent[k][k] = sqrt(sum);
  787. -
  788. -        /* set values of column k */
  789. -        for ( i=k+1; i<n; i++ )
  790. -        {
  791. -            sum = A_ent[i][k];
  792. -            A_piv = A_ent[k];
  793. -            A_row = A_ent[i];
  794. -            sum -= __ip__(A_row,A_piv,(int)k);
  795. -            /************************************************
  796. -            for ( j=0; j<k; j++ )
  797. -                sum -= A_ent[i][j]*A_ent[k][j];
  798. -                sum -= (*A_row++)*(*A_piv++);
  799. -            ************************************************/
  800. -            A_ent[j][i] = A_ent[i][j] = sum/A_ent[k][k];
  801. -        }
  802. -    }
  803. -
  804. -    return (A);
  805. -}
  806. //GO.SYSIN DD chfactor.c
  807. echo qrfactor.c 1>&2
  808. sed >qrfactor.c <<'//GO.SYSIN DD qrfactor.c' 's/^-//'
  809. -
  810. -/**************************************************************************
  811. -**
  812. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  813. -**
  814. -**                 Meschach Library
  815. -** 
  816. -** This Meschach Library is provided "as is" without any express 
  817. -** or implied warranty of any kind with respect to this software. 
  818. -** In particular the authors shall not be liable for any direct, 
  819. -** indirect, special, incidental or consequential damages arising 
  820. -** in any way from use of the software.
  821. -** 
  822. -** Everyone is granted permission to copy, modify and redistribute this
  823. -** Meschach Library, provided:
  824. -**  1.  All copies contain this copyright notice.
  825. -**  2.  All modified copies shall carry a notice stating who
  826. -**      made the last modification and the date of such modification.
  827. -**  3.  No charge is made for this software or works derived from it.  
  828. -**      This clause shall not be construed as constraining other software
  829. -**      distributed on the same medium as this software, nor is a
  830. -**      distribution fee considered a charge.
  831. -**
  832. -***************************************************************************/
  833. -
  834. -
  835. -/*
  836. -  This file contains the routines needed to perform QR factorisation
  837. -  of matrices, as well as Householder transformations.
  838. -  The internal "factored form" of a matrix A is not quite standard.
  839. -  The diagonal of A is replaced by the diagonal of R -- not by the 1st non-zero
  840. -  entries of the Householder vectors. The 1st non-zero entries are held in
  841. -  the diag parameter of QRfactor(). The reason for this non-standard
  842. -  representation is that it enables direct use of the Usolve() function
  843. -  rather than requiring that  a seperate function be written just for this case.
  844. -  See, e.g., QRsolve() below for more details.
  845. -  
  846. -*/
  847. -
  848. -
  849. -static    char    rcsid[] = "$Id: qrfactor.c,v 1.5 1994/01/13 05:35:07 des Exp $";
  850. -
  851. -#include    <stdio.h>
  852. -#include    <math.h>
  853. -#include        "matrix2.h"
  854. -
  855. -
  856. -
  857. -
  858. -
  859. -#define        sign(x)    ((x) > 0.0 ? 1 : ((x) < 0.0 ? -1 : 0 ))
  860. -
  861. -extern    VEC    *Usolve();    /* See matrix2.h */
  862. -
  863. -/* Note: The usual representation of a Householder transformation is taken
  864. -   to be:
  865. -   P = I - beta.u.uT
  866. -   where beta = 2/(uT.u) and u is called the Householder vector
  867. -   */
  868. -
  869. -/* QRfactor -- forms the QR factorisation of A -- factorisation stored in
  870. -   compact form as described above ( not quite standard format ) */
  871. -/* MAT    *QRfactor(A,diag,beta) */
  872. -MAT    *QRfactor(A,diag)
  873. -MAT    *A;
  874. -VEC    *diag /* ,*beta */;
  875. -{
  876. -    u_int    k,limit;
  877. -    Real    beta;
  878. -    static    VEC    *tmp1=VNULL;
  879. -    
  880. -    if ( ! A || ! diag )
  881. -    error(E_NULL,"QRfactor");
  882. -    limit = min(A->m,A->n);
  883. -    if ( diag->dim < limit )
  884. -    error(E_SIZES,"QRfactor");
  885. -    
  886. -    tmp1 = v_resize(tmp1,A->m);
  887. -    MEM_STAT_REG(tmp1,TYPE_VEC);
  888. -    
  889. -    for ( k=0; k<limit; k++ )
  890. -    {
  891. -    /* get H/holder vector for the k-th column */
  892. -    get_col(A,k,tmp1);
  893. -    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
  894. -    hhvec(tmp1,k,&beta,tmp1,&A->me[k][k]);
  895. -    diag->ve[k] = tmp1->ve[k];
  896. -    
  897. -    /* apply H/holder vector to remaining columns */
  898. -    /* hhtrcols(A,k,k+1,tmp1,beta->ve[k]); */
  899. -    hhtrcols(A,k,k+1,tmp1,beta);
  900. -    }
  901. -
  902. -    return (A);
  903. -}
  904. -
  905. -/* QRCPfactor -- forms the QR factorisation of A with column pivoting
  906. -   -- factorisation stored in compact form as described above
  907. -   ( not quite standard format )                */
  908. -/* MAT    *QRCPfactor(A,diag,beta,px) */
  909. -MAT    *QRCPfactor(A,diag,px)
  910. -MAT    *A;
  911. -VEC    *diag /* , *beta */;
  912. -PERM    *px;
  913. -{
  914. -    u_int    i, i_max, j, k, limit;
  915. -    static    VEC    *gamma=VNULL, *tmp1=VNULL, *tmp2=VNULL;
  916. -    Real    beta, maxgamma, sum, tmp;
  917. -    
  918. -    if ( ! A || ! diag || ! px )
  919. -    error(E_NULL,"QRCPfactor");
  920. -    limit = min(A->m,A->n);
  921. -    if ( diag->dim < limit || px->size != A->n )
  922. -    error(E_SIZES,"QRCPfactor");
  923. -    
  924. -    tmp1 = v_resize(tmp1,A->m);
  925. -    tmp2 = v_resize(tmp2,A->m);
  926. -    gamma = v_resize(gamma,A->n);
  927. -    MEM_STAT_REG(tmp1,TYPE_VEC);
  928. -    MEM_STAT_REG(tmp2,TYPE_VEC);
  929. -    MEM_STAT_REG(gamma,TYPE_VEC);
  930. -    
  931. -    /* initialise gamma and px */
  932. -    for ( j=0; j<A->n; j++ )
  933. -    {
  934. -    px->pe[j] = j;
  935. -    sum = 0.0;
  936. -    for ( i=0; i<A->m; i++ )
  937. -        sum += square(A->me[i][j]);
  938. -    gamma->ve[j] = sum;
  939. -    }
  940. -    
  941. -    for ( k=0; k<limit; k++ )
  942. -    {
  943. -    /* find "best" column to use */
  944. -    i_max = k;    maxgamma = gamma->ve[k];
  945. -    for ( i=k+1; i<A->n; i++ )
  946. -        /* Loop invariant:maxgamma=gamma[i_max]
  947. -           >=gamma[l];l=k,...,i-1 */
  948. -        if ( gamma->ve[i] > maxgamma )
  949. -        {    maxgamma = gamma->ve[i]; i_max = i;    }
  950. -    
  951. -    /* swap columns if necessary */
  952. -    if ( i_max != k )
  953. -    {
  954. -        /* swap gamma values */
  955. -        tmp = gamma->ve[k];
  956. -        gamma->ve[k] = gamma->ve[i_max];
  957. -        gamma->ve[i_max] = tmp;
  958. -        
  959. -        /* update column permutation */
  960. -        px_transp(px,k,i_max);
  961. -        
  962. -        /* swap columns of A */
  963. -        for ( i=0; i<A->m; i++ )
  964. -        {
  965. -        tmp = A->me[i][k];
  966. -        A->me[i][k] = A->me[i][i_max];
  967. -        A->me[i][i_max] = tmp;
  968. -        }
  969. -    }
  970. -    
  971. -    /* get H/holder vector for the k-th column */
  972. -    get_col(A,k,tmp1);
  973. -    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
  974. -    hhvec(tmp1,k,&beta,tmp1,&A->me[k][k]);
  975. -    diag->ve[k] = tmp1->ve[k];
  976. -    
  977. -    /* apply H/holder vector to remaining columns */
  978. -    /* hhtrcols(A,k,k+1,tmp1,beta->ve[k]); */
  979. -    hhtrcols(A,k,k+1,tmp1,beta);
  980. -    
  981. -    /* update gamma values */
  982. -    for ( j=k+1; j<A->n; j++ )
  983. -        gamma->ve[j] -= square(A->me[k][j]);
  984. -    }
  985. -
  986. -    return (A);
  987. -}
  988. -
  989. -/* Qsolve -- solves Qx = b, Q is an orthogonal matrix stored in compact
  990. -   form a la QRfactor() -- may be in-situ */
  991. -/* VEC    *_Qsolve(QR,diag,beta,b,x,tmp) */
  992. -VEC    *_Qsolve(QR,diag,b,x,tmp)
  993. -MAT    *QR;
  994. -VEC    *diag /* ,*beta */ , *b, *x, *tmp;
  995. -{
  996. -    u_int    dynamic;
  997. -    int        k, limit;
  998. -    Real    beta, r_ii, tmp_val;
  999. -    
  1000. -    limit = min(QR->m,QR->n);
  1001. -    dynamic = FALSE;
  1002. -    if ( ! QR || ! diag || ! b )
  1003. -    error(E_NULL,"_Qsolve");
  1004. -    if ( diag->dim < limit || b->dim != QR->m )
  1005. -    error(E_SIZES,"_Qsolve");
  1006. -    x = v_resize(x,QR->m);
  1007. -    if ( tmp == VNULL )
  1008. -    dynamic = TRUE;
  1009. -    tmp = v_resize(tmp,QR->m);
  1010. -    
  1011. -    /* apply H/holder transforms in normal order */
  1012. -    x = v_copy(b,x);
  1013. -    for ( k = 0 ; k < limit ; k++ )
  1014. -    {
  1015. -    get_col(QR,k,tmp);
  1016. -    r_ii = fabs(tmp->ve[k]);
  1017. -    tmp->ve[k] = diag->ve[k];
  1018. -    tmp_val = (r_ii*fabs(diag->ve[k]));
  1019. -    beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  1020. -    /* hhtrvec(tmp,beta->ve[k],k,x,x); */
  1021. -    hhtrvec(tmp,beta,k,x,x);
  1022. -    }
  1023. -    
  1024. -    if ( dynamic )
  1025. -    V_FREE(tmp);
  1026. -    
  1027. -    return (x);
  1028. -}
  1029. -
  1030. -/* makeQ -- constructs orthogonal matrix from Householder vectors stored in
  1031. -   compact QR form */
  1032. -/* MAT    *makeQ(QR,diag,beta,Qout) */
  1033. -MAT    *makeQ(QR,diag,Qout)
  1034. -MAT    *QR,*Qout;
  1035. -VEC    *diag /* , *beta */;
  1036. -{
  1037. -    static    VEC    *tmp1=VNULL,*tmp2=VNULL;
  1038. -    u_int    i, limit;
  1039. -    Real    beta, r_ii, tmp_val;
  1040. -    int    j;
  1041. -    
  1042. -    limit = min(QR->m,QR->n);
  1043. -    if ( ! QR || ! diag )
  1044. -    error(E_NULL,"makeQ");
  1045. -    if ( diag->dim < limit )
  1046. -    error(E_SIZES,"makeQ");
  1047. -    if ( Qout==(MAT *)NULL || Qout->m < QR->m || Qout->n < QR->m )
  1048. -    Qout = m_get(QR->m,QR->m);
  1049. -    
  1050. -    tmp1 = v_resize(tmp1,QR->m);    /* contains basis vec & columns of Q */
  1051. -    tmp2 = v_resize(tmp2,QR->m);    /* contains H/holder vectors */
  1052. -    MEM_STAT_REG(tmp1,TYPE_VEC);
  1053. -    MEM_STAT_REG(tmp2,TYPE_VEC);
  1054. -    
  1055. -    for ( i=0; i<QR->m ; i++ )
  1056. -    {    /* get i-th column of Q */
  1057. -    /* set up tmp1 as i-th basis vector */
  1058. -    for ( j=0; j<QR->m ; j++ )
  1059. -        tmp1->ve[j] = 0.0;
  1060. -    tmp1->ve[i] = 1.0;
  1061. -    
  1062. -    /* apply H/h transforms in reverse order */
  1063. -    for ( j=limit-1; j>=0; j-- )
  1064. -    {
  1065. -        get_col(QR,j,tmp2);
  1066. -        r_ii = fabs(tmp2->ve[j]);
  1067. -        tmp2->ve[j] = diag->ve[j];
  1068. -        tmp_val = (r_ii*fabs(diag->ve[j]));
  1069. -        beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  1070. -        /* hhtrvec(tmp2,beta->ve[j],j,tmp1,tmp1); */
  1071. -        hhtrvec(tmp2,beta,j,tmp1,tmp1);
  1072. -    }
  1073. -    
  1074. -    /* insert into Q */
  1075. -    set_col(Qout,i,tmp1);
  1076. -    }
  1077. -
  1078. -    return (Qout);
  1079. -}
  1080. -
  1081. -/* makeR -- constructs upper triangular matrix from QR (compact form)
  1082. -   -- may be in-situ (all it does is zero the lower 1/2) */
  1083. -MAT    *makeR(QR,Rout)
  1084. -MAT    *QR,*Rout;
  1085. -{
  1086. -    u_int    i,j;
  1087. -    
  1088. -    if ( QR==(MAT *)NULL )
  1089. -    error(E_NULL,"makeR");
  1090. -    Rout = m_copy(QR,Rout);
  1091. -    
  1092. -    for ( i=1; i<QR->m; i++ )
  1093. -    for ( j=0; j<QR->n && j<i; j++ )
  1094. -        Rout->me[i][j] = 0.0;
  1095. -    
  1096. -    return (Rout);
  1097. -}
  1098. -
  1099. -/* QRsolve -- solves the system Q.R.x=b where Q & R are stored in compact form
  1100. -   -- returns x, which is created if necessary */
  1101. -/* VEC    *QRsolve(QR,diag,beta,b,x) */
  1102. -VEC    *QRsolve(QR,diag,b,x)
  1103. -MAT    *QR;
  1104. -VEC    *diag /* , *beta */ , *b, *x;
  1105. -{
  1106. -    int    limit;
  1107. -    static    VEC    *tmp = VNULL;
  1108. -    
  1109. -    if ( ! QR || ! diag || ! b )
  1110. -    error(E_NULL,"QRsolve");
  1111. -    limit = min(QR->m,QR->n);
  1112. -    if ( diag->dim < limit || b->dim != QR->m )
  1113. -    error(E_SIZES,"QRsolve");
  1114. -    tmp = v_resize(tmp,limit);
  1115. -    MEM_STAT_REG(tmp,TYPE_VEC);
  1116. -
  1117. -    x = v_resize(x,QR->n);
  1118. -    _Qsolve(QR,diag,b,x,tmp);
  1119. -    x = Usolve(QR,x,x,0.0);
  1120. -    v_resize(x,QR->n);
  1121. -
  1122. -    return x;
  1123. -}
  1124. -
  1125. -/* QRCPsolve -- solves A.x = b where A is factored by QRCPfactor()
  1126. -   -- assumes that A is in the compact factored form */
  1127. -/* VEC    *QRCPsolve(QR,diag,beta,pivot,b,x) */
  1128. -VEC    *QRCPsolve(QR,diag,pivot,b,x)
  1129. -MAT    *QR;
  1130. -VEC    *diag /* , *beta */;
  1131. -PERM    *pivot;
  1132. -VEC    *b, *x;
  1133. -{
  1134. -    static    VEC    *tmp=VNULL;
  1135. -    
  1136. -    if ( ! QR || ! diag || ! pivot || ! b )
  1137. -    error(E_NULL,"QRCPsolve");
  1138. -    if ( (QR->m > diag->dim &&QR->n > diag->dim) || QR->n != pivot->size )
  1139. -    error(E_SIZES,"QRCPsolve");
  1140. -    
  1141. -    tmp = QRsolve(QR,diag /* , beta */ ,b,tmp);
  1142. -    MEM_STAT_REG(tmp,TYPE_VEC);
  1143. -    x = pxinv_vec(pivot,tmp,x);
  1144. -
  1145. -    return x;
  1146. -}
  1147. -
  1148. -/* Umlt -- compute out = upper_triang(U).x
  1149. -    -- may be in situ */
  1150. -static    VEC    *Umlt(U,x,out)
  1151. -MAT    *U;
  1152. -VEC    *x, *out;
  1153. -{
  1154. -    int        i, limit;
  1155. -
  1156. -    if ( U == MNULL || x == VNULL )
  1157. -    error(E_NULL,"Umlt");
  1158. -    limit = min(U->m,U->n);
  1159. -    if ( limit != x->dim )
  1160. -    error(E_SIZES,"Umlt");
  1161. -    if ( out == VNULL || out->dim < limit )
  1162. -    out = v_resize(out,limit);
  1163. -
  1164. -    for ( i = 0; i < limit; i++ )
  1165. -    out->ve[i] = __ip__(&(x->ve[i]),&(U->me[i][i]),limit - i);
  1166. -    return out;
  1167. -}
  1168. -
  1169. -/* UTmlt -- returns out = upper_triang(U)^T.x */
  1170. -static    VEC    *UTmlt(U,x,out)
  1171. -MAT    *U;
  1172. -VEC    *x, *out;
  1173. -{
  1174. -    Real    sum;
  1175. -    int        i, j, limit;
  1176. -
  1177. -    if ( U == MNULL || x == VNULL )
  1178. -    error(E_NULL,"UTmlt");
  1179. -    limit = min(U->m,U->n);
  1180. -    if ( out == VNULL || out->dim < limit )
  1181. -    out = v_resize(out,limit);
  1182. -
  1183. -    for ( i = limit-1; i >= 0; i-- )
  1184. -    {
  1185. -    sum = 0.0;
  1186. -    for ( j = 0; j <= i; j++ )
  1187. -        sum += U->me[j][i]*x->ve[j];
  1188. -    out->ve[i] = sum;
  1189. -    }
  1190. -    return out;
  1191. -}
  1192. -
  1193. -/* QRTsolve -- solve A^T.sc = c where the QR factors of A are stored in
  1194. -    compact form
  1195. -    -- returns sc
  1196. -    -- original due to Mike Osborne modified Wed 09th Dec 1992 */
  1197. -VEC *QRTsolve(A,diag,c,sc)
  1198. -MAT *A;
  1199. -VEC *diag, *c, *sc;
  1200. -{
  1201. -    int        i, j, k, n, p;
  1202. -    Real    beta, r_ii, s, tmp_val;
  1203. -
  1204. -    if ( ! A || ! diag || ! c )
  1205. -    error(E_NULL,"QRTsolve");
  1206. -    if ( diag->dim < min(A->m,A->n) )
  1207. -    error(E_SIZES,"QRTsolve");
  1208. -    sc = v_resize(sc,A->m);
  1209. -    n = sc->dim;
  1210. -    p = c->dim;
  1211. -    if ( n == p )
  1212. -    k = p-2;
  1213. -    else
  1214. -    k = p-1;
  1215. -    v_zero(sc);
  1216. -    sc->ve[0] = c->ve[0]/A->me[0][0];
  1217. -    if ( n ==  1)
  1218. -    return sc;
  1219. -    if ( p > 1)
  1220. -    {
  1221. -    for ( i = 1; i < p; i++ )
  1222. -    {
  1223. -        s = 0.0;
  1224. -        for ( j = 0; j < i; j++ )
  1225. -        s += A->me[j][i]*sc->ve[j];
  1226. -        if ( A->me[i][i] == 0.0 )
  1227. -        error(E_SING,"QRTsolve");
  1228. -        sc->ve[i]=(c->ve[i]-s)/A->me[i][i];
  1229. -    }
  1230. -    }
  1231. -    for (i = k; i >= 0; i--)
  1232. -    {
  1233. -    s = diag->ve[i]*sc->ve[i];
  1234. -    for ( j = i+1; j < n; j++ )
  1235. -        s += A->me[j][i]*sc->ve[j];
  1236. -    r_ii = fabs(A->me[i][i]);
  1237. -    tmp_val = (r_ii*fabs(diag->ve[i]));
  1238. -    beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  1239. -    tmp_val = beta*s;
  1240. -    sc->ve[i] -= tmp_val*diag->ve[i];
  1241. -    for ( j = i+1; j < n; j++ )
  1242. -        sc->ve[j] -= tmp_val*A->me[j][i];
  1243. -    }
  1244. -
  1245. -    return sc;
  1246. -}
  1247. -
  1248. -/* QRcondest -- returns an estimate of the 2-norm condition number of the
  1249. -        matrix factorised by QRfactor() or QRCPfactor()
  1250. -    -- note that as Q does not affect the 2-norm condition number,
  1251. -        it is not necessary to pass the diag, beta (or pivot) vectors
  1252. -    -- generates a lower bound on the true condition number
  1253. -    -- if the matrix is exactly singular, HUGE is returned
  1254. -    -- note that QRcondest() is likely to be more reliable for
  1255. -        matrices factored using QRCPfactor() */
  1256. -double    QRcondest(QR)
  1257. -MAT    *QR;
  1258. -{
  1259. -    static    VEC    *y=VNULL;
  1260. -    Real    norm1, norm2, sum, tmp1, tmp2;
  1261. -    int        i, j, limit;
  1262. -
  1263. -    if ( QR == MNULL )
  1264. -    error(E_NULL,"QRcondest");
  1265. -
  1266. -    limit = min(QR->m,QR->n);
  1267. -    for ( i = 0; i < limit; i++ )
  1268. -    if ( QR->me[i][i] == 0.0 )
  1269. -        return HUGE;
  1270. -
  1271. -    y = v_resize(y,limit);
  1272. -    MEM_STAT_REG(y,TYPE_VEC);
  1273. -    /* use the trick for getting a unit vector y with ||R.y||_inf small
  1274. -       from the LU condition estimator */
  1275. -    for ( i = 0; i < limit; i++ )
  1276. -    {
  1277. -    sum = 0.0;
  1278. -    for ( j = 0; j < i; j++ )
  1279. -        sum -= QR->me[j][i]*y->ve[j];
  1280. -    sum -= (sum < 0.0) ? 1.0 : -1.0;
  1281. -    y->ve[i] = sum / QR->me[i][i];
  1282. -    }
  1283. -    UTmlt(QR,y,y);
  1284. -
  1285. -    /* now apply inverse power method to R^T.R */
  1286. -    for ( i = 0; i < 3; i++ )
  1287. -    {
  1288. -    tmp1 = v_norm2(y);
  1289. -    sv_mlt(1/tmp1,y,y);
  1290. -    UTsolve(QR,y,y,0.0);
  1291. -    tmp2 = v_norm2(y);
  1292. -    sv_mlt(1/v_norm2(y),y,y);
  1293. -    Usolve(QR,y,y,0.0);
  1294. -    }
  1295. -    /* now compute approximation for ||R^{-1}||_2 */
  1296. -    norm1 = sqrt(tmp1)*sqrt(tmp2);
  1297. -
  1298. -    /* now use complementary approach to compute approximation to ||R||_2 */
  1299. -    for ( i = limit-1; i >= 0; i-- )
  1300. -    {
  1301. -    sum = 0.0;
  1302. -    for ( j = i+1; j < limit; j++ )
  1303. -        sum += QR->me[i][j]*y->ve[j];
  1304. -    y->ve[i] = (sum >= 0.0) ? 1.0 : -1.0;
  1305. -    y->ve[i] = (QR->me[i][i] >= 0.0) ? y->ve[i] : - y->ve[i];
  1306. -    }
  1307. -
  1308. -    /* now apply power method to R^T.R */
  1309. -    for ( i = 0; i < 3; i++ )
  1310. -    {
  1311. -    tmp1 = v_norm2(y);
  1312. -    sv_mlt(1/tmp1,y,y);
  1313. -    Umlt(QR,y,y);
  1314. -    tmp2 = v_norm2(y);
  1315. -    sv_mlt(1/tmp2,y,y);
  1316. -    UTmlt(QR,y,y);
  1317. -    }
  1318. -    norm2 = sqrt(tmp1)*sqrt(tmp2);
  1319. -
  1320. -    /* printf("QRcondest: norm1 = %g, norm2 = %g\n",norm1,norm2); */
  1321. -
  1322. -    return norm1*norm2;
  1323. -}
  1324. //GO.SYSIN DD qrfactor.c
  1325. echo solve.c 1>&2
  1326. sed >solve.c <<'//GO.SYSIN DD solve.c' 's/^-//'
  1327. -
  1328. -/**************************************************************************
  1329. -**
  1330. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  1331. -**
  1332. -**                 Meschach Library
  1333. -** 
  1334. -** This Meschach Library is provided "as is" without any express 
  1335. -** or implied warranty of any kind with respect to this software. 
  1336. -** In particular the authors shall not be liable for any direct, 
  1337. -** indirect, special, incidental or consequential damages arising 
  1338. -** in any way from use of the software.
  1339. -** 
  1340. -** Everyone is granted permission to copy, modify and redistribute this
  1341. -** Meschach Library, provided:
  1342. -**  1.  All copies contain this copyright notice.
  1343. -**  2.  All modified copies shall carry a notice stating who
  1344. -**      made the last modification and the date of such modification.
  1345. -**  3.  No charge is made for this software or works derived from it.  
  1346. -**      This clause shall not be construed as constraining other software
  1347. -**      distributed on the same medium as this software, nor is a
  1348. -**      distribution fee considered a charge.
  1349. -**
  1350. -***************************************************************************/
  1351. -
  1352. -
  1353. -/*
  1354. -    Matrix factorisation routines to work with the other matrix files.
  1355. -*/
  1356. -
  1357. -/* solve.c 1.2 11/25/87 */
  1358. -static    char    rcsid[] = "$Id: solve.c,v 1.3 1994/01/13 05:29:57 des Exp $";
  1359. -
  1360. -#include    <stdio.h>
  1361. -#include    <math.h>
  1362. -#include        "matrix2.h"
  1363. -
  1364. -
  1365. -
  1366. -
  1367. -
  1368. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  1369. -
  1370. -/* Usolve -- back substitution with optional over-riding diagonal
  1371. -        -- can be in-situ but doesn't need to be */
  1372. -VEC    *Usolve(matrix,b,out,diag)
  1373. -MAT    *matrix;
  1374. -VEC    *b, *out;
  1375. -double    diag;
  1376. -{
  1377. -    u_int    dim /* , j */;
  1378. -    int    i, i_lim;
  1379. -    Real    **mat_ent, *mat_row, *b_ent, *out_ent, *out_col, sum;
  1380. -
  1381. -    if ( matrix==(MAT *)NULL || b==(VEC *)NULL )
  1382. -        error(E_NULL,"Usolve");
  1383. -    dim = min(matrix->m,matrix->n);
  1384. -    if ( b->dim < dim )
  1385. -        error(E_SIZES,"Usolve");
  1386. -    if ( out==(VEC *)NULL || out->dim < dim )
  1387. -        out = v_resize(out,matrix->n);
  1388. -    mat_ent = matrix->me;    b_ent = b->ve;    out_ent = out->ve;
  1389. -
  1390. -    for ( i=dim-1; i>=0; i-- )
  1391. -        if ( b_ent[i] != 0.0 )
  1392. -            break;
  1393. -        else
  1394. -            out_ent[i] = 0.0;
  1395. -    i_lim = i;
  1396. -
  1397. -    for (    ; i>=0; i-- )
  1398. -    {
  1399. -        sum = b_ent[i];
  1400. -        mat_row = &(mat_ent[i][i+1]);
  1401. -        out_col = &(out_ent[i+1]);
  1402. -        sum -= __ip__(mat_row,out_col,i_lim-i);
  1403. -        /******************************************************
  1404. -        for ( j=i+1; j<=i_lim; j++ )
  1405. -            sum -= mat_ent[i][j]*out_ent[j];
  1406. -            sum -= (*mat_row++)*(*out_col++);
  1407. -        ******************************************************/
  1408. -        if ( diag==0.0 )
  1409. -        {
  1410. -            if ( mat_ent[i][i]==0.0 )
  1411. -                error(E_SING,"Usolve");
  1412. -            else
  1413. -                out_ent[i] = sum/mat_ent[i][i];
  1414. -        }
  1415. -        else
  1416. -            out_ent[i] = sum/diag;
  1417. -    }
  1418. -
  1419. -    return (out);
  1420. -}
  1421. -
  1422. -/* Lsolve -- forward elimination with (optional) default diagonal value */
  1423. -VEC    *Lsolve(matrix,b,out,diag)
  1424. -MAT    *matrix;
  1425. -VEC    *b,*out;
  1426. -double    diag;
  1427. -{
  1428. -    u_int    dim, i, i_lim /* , j */;
  1429. -    Real    **mat_ent, *mat_row, *b_ent, *out_ent, *out_col, sum;
  1430. -
  1431. -    if ( matrix==(MAT *)NULL || b==(VEC *)NULL )
  1432. -        error(E_NULL,"Lsolve");
  1433. -    dim = min(matrix->m,matrix->n);
  1434. -    if ( b->dim < dim )
  1435. -        error(E_SIZES,"Lsolve");
  1436. -    if ( out==(VEC *)NULL || out->dim < dim )
  1437. -        out = v_resize(out,matrix->n);
  1438. -    mat_ent = matrix->me;    b_ent = b->ve;    out_ent = out->ve;
  1439. -
  1440. -    for ( i=0; i<dim; i++ )
  1441. -        if ( b_ent[i] != 0.0 )
  1442. -            break;
  1443. -        else
  1444. -            out_ent[i] = 0.0;
  1445. -    i_lim = i;
  1446. -
  1447. -    for (    ; i<dim; i++ )
  1448. -    {
  1449. -        sum = b_ent[i];
  1450. -        mat_row = &(mat_ent[i][i_lim]);
  1451. -        out_col = &(out_ent[i_lim]);
  1452. -        sum -= __ip__(mat_row,out_col,(int)(i-i_lim));
  1453. -        /*****************************************************
  1454. -        for ( j=i_lim; j<i; j++ )
  1455. -            sum -= mat_ent[i][j]*out_ent[j];
  1456. -            sum -= (*mat_row++)*(*out_col++);
  1457. -        ******************************************************/
  1458. -        if ( diag==0.0 )
  1459. -        {
  1460. -            if ( mat_ent[i][i]==0.0 )
  1461. -                error(E_SING,"Lsolve");
  1462. -            else
  1463. -                out_ent[i] = sum/mat_ent[i][i];
  1464. -        }
  1465. -        else
  1466. -            out_ent[i] = sum/diag;
  1467. -    }
  1468. -
  1469. -    return (out);
  1470. -}
  1471. -
  1472. -
  1473. -/* UTsolve -- forward elimination with (optional) default diagonal value
  1474. -        using UPPER triangular part of matrix */
  1475. -VEC    *UTsolve(U,b,out,diag)
  1476. -MAT    *U;
  1477. -VEC    *b,*out;
  1478. -double    diag;
  1479. -{
  1480. -    u_int    dim, i, i_lim;
  1481. -    Real    **U_me, *b_ve, *out_ve, tmp, invdiag;
  1482. -    
  1483. -    if ( ! U || ! b )
  1484. -    error(E_NULL,"UTsolve");
  1485. -    dim = min(U->m,U->n);
  1486. -    if ( b->dim < dim )
  1487. -    error(E_SIZES,"UTsolve");
  1488. -    out = v_resize(out,U->n);
  1489. -    U_me = U->me;    b_ve = b->ve;    out_ve = out->ve;
  1490. -    
  1491. -    for ( i=0; i<dim; i++ )
  1492. -    if ( b_ve[i] != 0.0 )
  1493. -        break;
  1494. -    else
  1495. -        out_ve[i] = 0.0;
  1496. -    i_lim = i;
  1497. -    if ( b != out )
  1498. -    {
  1499. -    __zero__(out_ve,out->dim);
  1500. -    MEM_COPY(&(b_ve[i_lim]),&(out_ve[i_lim]),(dim-i_lim)*sizeof(Real));
  1501. -    }
  1502. -
  1503. -    if ( diag == 0.0 )
  1504. -    {
  1505. -    for (    ; i<dim; i++ )
  1506. -    {
  1507. -        tmp = U_me[i][i];
  1508. -        if ( tmp == 0.0 )
  1509. -        error(E_SING,"UTsolve");
  1510. -        out_ve[i] /= tmp;
  1511. -        __mltadd__(&(out_ve[i+1]),&(U_me[i][i+1]),-out_ve[i],dim-i-1);
  1512. -    }
  1513. -    }
  1514. -    else
  1515. -    {
  1516. -    invdiag = 1.0/diag;
  1517. -    for (    ; i<dim; i++ )
  1518. -    {
  1519. -        out_ve[i] *= invdiag;
  1520. -        __mltadd__(&(out_ve[i+1]),&(U_me[i][i+1]),-out_ve[i],dim-i-1);
  1521. -    }
  1522. -    }
  1523. -    return (out);
  1524. -}
  1525. -
  1526. -/* Dsolve -- solves Dx=b where D is the diagonal of A -- may be in-situ */
  1527. -VEC    *Dsolve(A,b,x)
  1528. -MAT    *A;
  1529. -VEC    *b,*x;
  1530. -{
  1531. -    u_int    dim, i;
  1532. -    
  1533. -    if ( ! A || ! b )
  1534. -    error(E_NULL,"Dsolve");
  1535. -    dim = min(A->m,A->n);
  1536. -    if ( b->dim < dim )
  1537. -    error(E_SIZES,"Dsolve");
  1538. -    x = v_resize(x,A->n);
  1539. -    
  1540. -    dim = b->dim;
  1541. -    for ( i=0; i<dim; i++ )
  1542. -    if ( A->me[i][i] == 0.0 )
  1543. -        error(E_SING,"Dsolve");
  1544. -    else
  1545. -        x->ve[i] = b->ve[i]/A->me[i][i];
  1546. -    
  1547. -    return (x);
  1548. -}
  1549. -
  1550. -/* LTsolve -- back substitution with optional over-riding diagonal
  1551. -        using the LOWER triangular part of matrix
  1552. -        -- can be in-situ but doesn't need to be */
  1553. -VEC    *LTsolve(L,b,out,diag)
  1554. -MAT    *L;
  1555. -VEC    *b, *out;
  1556. -double    diag;
  1557. -{
  1558. -    u_int    dim;
  1559. -    int        i, i_lim;
  1560. -    Real    **L_me, *b_ve, *out_ve, tmp, invdiag;
  1561. -    
  1562. -    if ( ! L || ! b )
  1563. -    error(E_NULL,"LTsolve");
  1564. -    dim = min(L->m,L->n);
  1565. -    if ( b->dim < dim )
  1566. -    error(E_SIZES,"LTsolve");
  1567. -    out = v_resize(out,L->n);
  1568. -    L_me = L->me;    b_ve = b->ve;    out_ve = out->ve;
  1569. -    
  1570. -    for ( i=dim-1; i>=0; i-- )
  1571. -    if ( b_ve[i] != 0.0 )
  1572. -        break;
  1573. -    i_lim = i;
  1574. -
  1575. -    if ( b != out )
  1576. -    {
  1577. -    __zero__(out_ve,out->dim);
  1578. -    MEM_COPY(b_ve,out_ve,(i_lim+1)*sizeof(Real));
  1579. -    }
  1580. -
  1581. -    if ( diag == 0.0 )
  1582. -    {
  1583. -    for (        ; i>=0; i-- )
  1584. -    {
  1585. -        tmp = L_me[i][i];
  1586. -        if ( tmp == 0.0 )
  1587. -        error(E_SING,"LTsolve");
  1588. -        out_ve[i] /= tmp;
  1589. -        __mltadd__(out_ve,L_me[i],-out_ve[i],i);
  1590. -    }
  1591. -    }
  1592. -    else
  1593. -    {
  1594. -    invdiag = 1.0/diag;
  1595. -    for (        ; i>=0; i-- )
  1596. -    {
  1597. -        out_ve[i] *= invdiag;
  1598. -        __mltadd__(out_ve,L_me[i],-out_ve[i],i);
  1599. -    }
  1600. -    }
  1601. -    
  1602. -    return (out);
  1603. -}
  1604. //GO.SYSIN DD solve.c
  1605. echo hsehldr.c 1>&2
  1606. sed >hsehldr.c <<'//GO.SYSIN DD hsehldr.c' 's/^-//'
  1607. -
  1608. -/**************************************************************************
  1609. -**
  1610. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  1611. -**
  1612. -**                 Meschach Library
  1613. -** 
  1614. -** This Meschach Library is provided "as is" without any express 
  1615. -** or implied warranty of any kind with respect to this software. 
  1616. -** In particular the authors shall not be liable for any direct, 
  1617. -** indirect, special, incidental or consequential damages arising 
  1618. -** in any way from use of the software.
  1619. -** 
  1620. -** Everyone is granted permission to copy, modify and redistribute this
  1621. -** Meschach Library, provided:
  1622. -**  1.  All copies contain this copyright notice.
  1623. -**  2.  All modified copies shall carry a notice stating who
  1624. -**      made the last modification and the date of such modification.
  1625. -**  3.  No charge is made for this software or works derived from it.  
  1626. -**      This clause shall not be construed as constraining other software
  1627. -**      distributed on the same medium as this software, nor is a
  1628. -**      distribution fee considered a charge.
  1629. -**
  1630. -***************************************************************************/
  1631. -
  1632. -
  1633. -/*
  1634. -        Files for matrix computations
  1635. -
  1636. -    Householder transformation file. Contains routines for calculating
  1637. -    householder transformations, applying them to vectors and matrices
  1638. -    by both row & column.
  1639. -*/
  1640. -
  1641. -/* hsehldr.c 1.3 10/8/87 */
  1642. -static    char    rcsid[] = "$Id: hsehldr.c,v 1.2 1994/01/13 05:36:29 des Exp $";
  1643. -
  1644. -#include    <stdio.h>
  1645. -#include    <math.h>
  1646. -#include    "matrix.h"
  1647. -#include        "matrix2.h"
  1648. -
  1649. -
  1650. -/* hhvec -- calulates Householder vector to eliminate all entries after the
  1651. -    i0 entry of the vector vec. It is returned as out. May be in-situ */
  1652. -VEC    *hhvec(vec,i0,beta,out,newval)
  1653. -VEC    *vec,*out;
  1654. -u_int    i0;
  1655. -Real    *beta,*newval;
  1656. -{
  1657. -    Real    norm;
  1658. -
  1659. -    out = _v_copy(vec,out,i0);
  1660. -    norm = sqrt(_in_prod(out,out,i0));
  1661. -    if ( norm <= 0.0 )
  1662. -    {
  1663. -        *beta = 0.0;
  1664. -        return (out);
  1665. -    }
  1666. -    *beta = 1.0/(norm * (norm+fabs(out->ve[i0])));
  1667. -    if ( out->ve[i0] > 0.0 )
  1668. -        *newval = -norm;
  1669. -    else
  1670. -        *newval = norm;
  1671. -    out->ve[i0] -= *newval;
  1672. -
  1673. -    return (out);
  1674. -}
  1675. -
  1676. -/* hhtrvec -- apply Householder transformation to vector -- may be in-situ */
  1677. -VEC    *hhtrvec(hh,beta,i0,in,out)
  1678. -VEC    *hh,*in,*out;    /* hh = Householder vector */
  1679. -u_int    i0;
  1680. -double    beta;
  1681. -{
  1682. -    Real    scale;
  1683. -    /* u_int    i; */
  1684. -
  1685. -    if ( hh==(VEC *)NULL || in==(VEC *)NULL )
  1686. -        error(E_NULL,"hhtrvec");
  1687. -    if ( in->dim != hh->dim )
  1688. -        error(E_SIZES,"hhtrvec");
  1689. -    if ( i0 > in->dim )
  1690. -        error(E_BOUNDS,"hhtrvec");
  1691. -
  1692. -    scale = beta*_in_prod(hh,in,i0);
  1693. -    out = v_copy(in,out);
  1694. -    __mltadd__(&(out->ve[i0]),&(hh->ve[i0]),-scale,(int)(in->dim-i0));
  1695. -    /************************************************************
  1696. -    for ( i=i0; i<in->dim; i++ )
  1697. -        out->ve[i] = in->ve[i] - scale*hh->ve[i];
  1698. -    ************************************************************/
  1699. -
  1700. -    return (out);
  1701. -}
  1702. -
  1703. -/* hhtrrows -- transform a matrix by a Householder vector by rows
  1704. -    starting at row i0 from column j0 -- in-situ */
  1705. -MAT    *hhtrrows(M,i0,j0,hh,beta)
  1706. -MAT    *M;
  1707. -u_int    i0, j0;
  1708. -VEC    *hh;
  1709. -double    beta;
  1710. -{
  1711. -    Real    ip, scale;
  1712. -    int    i /*, j */;
  1713. -
  1714. -    if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  1715. -        error(E_NULL,"hhtrrows");
  1716. -    if ( M->n != hh->dim )
  1717. -        error(E_RANGE,"hhtrrows");
  1718. -    if ( i0 > M->m || j0 > M->n )
  1719. -        error(E_BOUNDS,"hhtrrows");
  1720. -
  1721. -    if ( beta == 0.0 )    return (M);
  1722. -
  1723. -    /* for each row ... */
  1724. -    for ( i = i0; i < M->m; i++ )
  1725. -    {    /* compute inner product */
  1726. -        ip = __ip__(&(M->me[i][j0]),&(hh->ve[j0]),(int)(M->n-j0));
  1727. -        /**************************************************
  1728. -        ip = 0.0;
  1729. -        for ( j = j0; j < M->n; j++ )
  1730. -            ip += M->me[i][j]*hh->ve[j];
  1731. -        **************************************************/
  1732. -        scale = beta*ip;
  1733. -        if ( scale == 0.0 )
  1734. -            continue;
  1735. -
  1736. -        /* do operation */
  1737. -        __mltadd__(&(M->me[i][j0]),&(hh->ve[j0]),-scale,
  1738. -                            (int)(M->n-j0));
  1739. -        /**************************************************
  1740. -        for ( j = j0; j < M->n; j++ )
  1741. -            M->me[i][j] -= scale*hh->ve[j];
  1742. -        **************************************************/
  1743. -    }
  1744. -
  1745. -    return (M);
  1746. -}
  1747. -
  1748. -
  1749. -/* hhtrcols -- transform a matrix by a Householder vector by columns
  1750. -    starting at row i0 from column j0 -- in-situ */
  1751. -MAT    *hhtrcols(M,i0,j0,hh,beta)
  1752. -MAT    *M;
  1753. -u_int    i0, j0;
  1754. -VEC    *hh;
  1755. -double    beta;
  1756. -{
  1757. -    /* Real    ip, scale; */
  1758. -    int    i /*, k */;
  1759. -    static    VEC    *w = VNULL;
  1760. -
  1761. -    if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  1762. -        error(E_NULL,"hhtrcols");
  1763. -    if ( M->m != hh->dim )
  1764. -        error(E_SIZES,"hhtrcols");
  1765. -    if ( i0 > M->m || j0 > M->n )
  1766. -        error(E_BOUNDS,"hhtrcols");
  1767. -
  1768. -    if ( beta == 0.0 )    return (M);
  1769. -
  1770. -    w = v_resize(w,M->n);
  1771. -    MEM_STAT_REG(w,TYPE_VEC);
  1772. -    v_zero(w);
  1773. -
  1774. -    for ( i = i0; i < M->m; i++ )
  1775. -        if ( hh->ve[i] != 0.0 )
  1776. -        __mltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  1777. -                            (int)(M->n-j0));
  1778. -    for ( i = i0; i < M->m; i++ )
  1779. -        if ( hh->ve[i] != 0.0 )
  1780. -        __mltadd__(&(M->me[i][j0]),&(w->ve[j0]),-beta*hh->ve[i],
  1781. -                            (int)(M->n-j0));
  1782. -    return (M);
  1783. -}
  1784. -
  1785. //GO.SYSIN DD hsehldr.c
  1786. echo givens.c 1>&2
  1787. sed >givens.c <<'//GO.SYSIN DD givens.c' 's/^-//'
  1788. -
  1789. -/**************************************************************************
  1790. -**
  1791. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  1792. -**
  1793. -**                 Meschach Library
  1794. -** 
  1795. -** This Meschach Library is provided "as is" without any express 
  1796. -** or implied warranty of any kind with respect to this software. 
  1797. -** In particular the authors shall not be liable for any direct, 
  1798. -** indirect, special, incidental or consequential damages arising 
  1799. -** in any way from use of the software.
  1800. -** 
  1801. -** Everyone is granted permission to copy, modify and redistribute this
  1802. -** Meschach Library, provided:
  1803. -**  1.  All copies contain this copyright notice.
  1804. -**  2.  All modified copies shall carry a notice stating who
  1805. -**      made the last modification and the date of such modification.
  1806. -**  3.  No charge is made for this software or works derived from it.  
  1807. -**      This clause shall not be construed as constraining other software
  1808. -**      distributed on the same medium as this software, nor is a
  1809. -**      distribution fee considered a charge.
  1810. -**
  1811. -***************************************************************************/
  1812. -
  1813. -
  1814. -
  1815. -/*
  1816. -        Files for matrix computations
  1817. -
  1818. -    Givens operations file. Contains routines for calculating and
  1819. -    applying givens rotations for/to vectors and also to matrices by
  1820. -    row and by column.
  1821. -*/
  1822. -
  1823. -/* givens.c 1.2 11/25/87 */
  1824. -static    char    rcsid[] = "$Id: givens.c,v 1.2 1994/01/13 05:39:42 des Exp $";
  1825. -
  1826. -#include    <stdio.h>
  1827. -#include    <math.h>
  1828. -#include    "matrix.h"
  1829. -#include        "matrix2.h"
  1830. -
  1831. -/* givens -- returns c,s parameters for Givens rotation to
  1832. -        eliminate y in the vector [ x y ]' */
  1833. -void    givens(x,y,c,s)
  1834. -double  x,y;
  1835. -Real    *c,*s;
  1836. -{
  1837. -    Real    norm;
  1838. -
  1839. -    norm = sqrt(x*x+y*y);
  1840. -    if ( norm == 0.0 )
  1841. -    {    *c = 1.0;    *s = 0.0;    }    /* identity */
  1842. -    else
  1843. -    {    *c = x/norm;    *s = y/norm;    }
  1844. -}
  1845. -
  1846. -/* rot_vec -- apply Givens rotation to x's i & k components */
  1847. -VEC    *rot_vec(x,i,k,c,s,out)
  1848. -VEC    *x,*out;
  1849. -u_int    i,k;
  1850. -double    c,s;
  1851. -{
  1852. -    Real    temp;
  1853. -
  1854. -    if ( x==VNULL )
  1855. -        error(E_NULL,"rot_vec");
  1856. -    if ( i >= x->dim || k >= x->dim )
  1857. -        error(E_RANGE,"rot_vec");
  1858. -    out = v_copy(x,out);
  1859. -
  1860. -    /* temp = c*out->ve[i] + s*out->ve[k]; */
  1861. -    temp = c*v_entry(out,i) + s*v_entry(out,k);
  1862. -    /* out->ve[k] = -s*out->ve[i] + c*out->ve[k]; */
  1863. -    v_set_val(out,k,-s*v_entry(out,i)+c*v_entry(out,k));
  1864. -    /* out->ve[i] = temp; */
  1865. -    v_set_val(out,i,temp);
  1866. -
  1867. -    return (out);
  1868. -}
  1869. -
  1870. -/* rot_rows -- premultiply mat by givens rotation described by c,s */
  1871. -MAT    *rot_rows(mat,i,k,c,s,out)
  1872. -MAT    *mat,*out;
  1873. -u_int    i,k;
  1874. -double    c,s;
  1875. -{
  1876. -    u_int    j;
  1877. -    Real    temp;
  1878. -
  1879. -    if ( mat==(MAT *)NULL )
  1880. -        error(E_NULL,"rot_rows");
  1881. -    if ( i >= mat->m || k >= mat->m )
  1882. -        error(E_RANGE,"rot_rows");
  1883. -    out = m_copy(mat,out);
  1884. -
  1885. -    for ( j=0; j<mat->n; j++ )
  1886. -    {
  1887. -        /* temp = c*out->me[i][j] + s*out->me[k][j]; */
  1888. -        temp = c*m_entry(out,i,j) + s*m_entry(out,k,j);
  1889. -        /* out->me[k][j] = -s*out->me[i][j] + c*out->me[k][j]; */
  1890. -        m_set_val(out,k,j, -s*m_entry(out,i,j) + c*m_entry(out,k,j));
  1891. -        /* out->me[i][j] = temp; */
  1892. -        m_set_val(out,i,j, temp);
  1893. -    }
  1894. -
  1895. -    return (out);
  1896. -}
  1897. -
  1898. -/* rot_cols -- postmultiply mat by givens rotation described by c,s */
  1899. -MAT    *rot_cols(mat,i,k,c,s,out)
  1900. -MAT    *mat,*out;
  1901. -u_int    i,k;
  1902. -double    c,s;
  1903. -{
  1904. -    u_int    j;
  1905. -    Real    temp;
  1906. -
  1907. -    if ( mat==(MAT *)NULL )
  1908. -        error(E_NULL,"rot_cols");
  1909. -    if ( i >= mat->n || k >= mat->n )
  1910. -        error(E_RANGE,"rot_cols");
  1911. -    out = m_copy(mat,out);
  1912. -
  1913. -    for ( j=0; j<mat->m; j++ )
  1914. -    {
  1915. -        /* temp = c*out->me[j][i] + s*out->me[j][k]; */
  1916. -        temp = c*m_entry(out,j,i) + s*m_entry(out,j,k);
  1917. -        /* out->me[j][k] = -s*out->me[j][i] + c*out->me[j][k]; */
  1918. -        m_set_val(out,j,k, -s*m_entry(out,j,i) + c*m_entry(out,j,k));
  1919. -        /* out->me[j][i] = temp; */
  1920. -        m_set_val(out,j,i,temp);
  1921. -    }
  1922. -
  1923. -    return (out);
  1924. -}
  1925. -
  1926. //GO.SYSIN DD givens.c
  1927. echo update.c 1>&2
  1928. sed >update.c <<'//GO.SYSIN DD update.c' 's/^-//'
  1929. -
  1930. -/**************************************************************************
  1931. -**
  1932. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  1933. -**
  1934. -**                 Meschach Library
  1935. -** 
  1936. -** This Meschach Library is provided "as is" without any express 
  1937. -** or implied warranty of any kind with respect to this software. 
  1938. -** In particular the authors shall not be liable for any direct, 
  1939. -** indirect, special, incidental or consequential damages arising 
  1940. -** in any way from use of the software.
  1941. -** 
  1942. -** Everyone is granted permission to copy, modify and redistribute this
  1943. -** Meschach Library, provided:
  1944. -**  1.  All copies contain this copyright notice.
  1945. -**  2.  All modified copies shall carry a notice stating who
  1946. -**      made the last modification and the date of such modification.
  1947. -**  3.  No charge is made for this software or works derived from it.  
  1948. -**      This clause shall not be construed as constraining other software
  1949. -**      distributed on the same medium as this software, nor is a
  1950. -**      distribution fee considered a charge.
  1951. -**
  1952. -***************************************************************************/
  1953. -
  1954. -
  1955. -/*
  1956. -    Matrix factorisation routines to work with the other matrix files.
  1957. -*/
  1958. -
  1959. -/* update.c 1.3 11/25/87 */
  1960. -static    char    rcsid[] = "$Id: update.c,v 1.2 1994/01/13 05:26:06 des Exp $";
  1961. -
  1962. -#include    <stdio.h>
  1963. -#include    <math.h>
  1964. -#include    "matrix.h"
  1965. -#include        "matrix2.h"
  1966. -
  1967. -
  1968. -
  1969. -
  1970. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  1971. -
  1972. -/* LDLupdate -- updates a CHolesky factorisation, replacing LDL' by
  1973. -    MD~M' = LDL' + alpha.w.w' Note: w is overwritten
  1974. -    Ref: Gill et al Math Comp 28, p516 Algorithm C1 */
  1975. -MAT    *LDLupdate(CHmat,w,alpha)
  1976. -MAT    *CHmat;
  1977. -VEC    *w;
  1978. -double    alpha;
  1979. -{
  1980. -    u_int    i,j;
  1981. -    Real    diag,new_diag,beta,p;
  1982. -
  1983. -    if ( CHmat==(MAT *)NULL || w==(VEC *)NULL )
  1984. -        error(E_NULL,"LDLupdate");
  1985. -    if ( CHmat->m != CHmat->n || w->dim != CHmat->m )
  1986. -        error(E_SIZES,"LDLupdate");
  1987. -
  1988. -    for ( j=0; j < w->dim; j++ )
  1989. -    {
  1990. -        p = w->ve[j];
  1991. -        diag = CHmat->me[j][j];
  1992. -        new_diag = CHmat->me[j][j] = diag + alpha*p*p;
  1993. -        if ( new_diag <= 0.0 )
  1994. -            error(E_POSDEF,"LDLupdate");
  1995. -        beta = p*alpha/new_diag;
  1996. -        alpha *= diag/new_diag;
  1997. -
  1998. -        for ( i=j+1; i < w->dim; i++ )
  1999. -        {
  2000. -            w->ve[i] -= p*CHmat->me[i][j];
  2001. -            CHmat->me[i][j] += beta*w->ve[i];
  2002. -            CHmat->me[j][i] = CHmat->me[i][j];
  2003. -        }
  2004. -    }
  2005. -
  2006. -    return (CHmat);
  2007. -}
  2008. -
  2009. -
  2010. -/* QRupdate -- updates QR factorisation in expanded form (seperate matrices)
  2011. -    Finds Q+, R+ s.t. Q+.R+ = Q.(R+u.v') and Q+ orthogonal, R+ upper triang
  2012. -    Ref: Golub & van Loan Matrix Computations pp437-443
  2013. -    -- does not update Q if it is NULL */
  2014. -MAT    *QRupdate(Q,R,u,v)
  2015. -MAT    *Q,*R;
  2016. -VEC    *u,*v;
  2017. -{
  2018. -    int    i,j,k;
  2019. -    Real    c,s,temp;
  2020. -
  2021. -    if ( ! R || ! u || ! v )
  2022. -        error(E_NULL,"QRupdate");
  2023. -    if ( ( Q && ( Q->m != Q->n || R->m != Q->n ) ) ||
  2024. -                    u->dim != R->m || v->dim != R->n )
  2025. -        error(E_SIZES,"QRupdate");
  2026. -
  2027. -    /* find largest k s.t. u[k] != 0 */
  2028. -    for ( k=R->m-1; k>=0; k-- )
  2029. -        if ( u->ve[k] != 0.0 )
  2030. -            break;
  2031. -
  2032. -    /* transform R+u.v' to Hessenberg form */
  2033. -    for ( i=k-1; i>=0; i-- )
  2034. -    {
  2035. -        /* get Givens rotation */
  2036. -        givens(u->ve[i],u->ve[i+1],&c,&s);
  2037. -        rot_rows(R,i,i+1,c,s,R);
  2038. -        if ( Q )
  2039. -            rot_cols(Q,i,i+1,c,s,Q);
  2040. -        rot_vec(u,i,i+1,c,s,u);
  2041. -    }
  2042. -
  2043. -    /* add into R */
  2044. -    temp = u->ve[0];
  2045. -    for ( j=0; j<R->n; j++ )
  2046. -        R->me[0][j] += temp*v->ve[j];
  2047. -
  2048. -    /* transform Hessenberg to upper triangular */
  2049. -    for ( i=0; i<k; i++ )
  2050. -    {
  2051. -        givens(R->me[i][i],R->me[i+1][i],&c,&s);
  2052. -        rot_rows(R,i,i+1,c,s,R);
  2053. -        if ( Q )
  2054. -            rot_cols(Q,i,i+1,c,s,Q);
  2055. -    }
  2056. -
  2057. -    return R;
  2058. -}
  2059. -
  2060. //GO.SYSIN DD update.c
  2061. echo norm.c 1>&2
  2062. sed >norm.c <<'//GO.SYSIN DD norm.c' 's/^-//'
  2063. -
  2064. -/**************************************************************************
  2065. -**
  2066. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2067. -**
  2068. -**                 Meschach Library
  2069. -** 
  2070. -** This Meschach Library is provided "as is" without any express 
  2071. -** or implied warranty of any kind with respect to this software. 
  2072. -** In particular the authors shall not be liable for any direct, 
  2073. -** indirect, special, incidental or consequential damages arising 
  2074. -** in any way from use of the software.
  2075. -** 
  2076. -** Everyone is granted permission to copy, modify and redistribute this
  2077. -** Meschach Library, provided:
  2078. -**  1.  All copies contain this copyright notice.
  2079. -**  2.  All modified copies shall carry a notice stating who
  2080. -**      made the last modification and the date of such modification.
  2081. -**  3.  No charge is made for this software or works derived from it.  
  2082. -**      This clause shall not be construed as constraining other software
  2083. -**      distributed on the same medium as this software, nor is a
  2084. -**      distribution fee considered a charge.
  2085. -**
  2086. -***************************************************************************/
  2087. -
  2088. -
  2089. -/*
  2090. -    A collection of functions for computing norms: scaled and unscaled
  2091. -*/
  2092. -static    char    rcsid[] = "$Id: norm.c,v 1.6 1994/01/13 05:34:35 des Exp $";
  2093. -
  2094. -#include    <stdio.h>
  2095. -#include    <math.h>
  2096. -#include    "matrix.h"
  2097. -
  2098. -
  2099. -/* _v_norm1 -- computes (scaled) 1-norms of vectors */
  2100. -double    _v_norm1(x,scale)
  2101. -VEC    *x, *scale;
  2102. -{
  2103. -    int    i, dim;
  2104. -    Real    s, sum;
  2105. -
  2106. -    if ( x == (VEC *)NULL )
  2107. -        error(E_NULL,"_v_norm1");
  2108. -    dim = x->dim;
  2109. -
  2110. -    sum = 0.0;
  2111. -    if ( scale == (VEC *)NULL )
  2112. -        for ( i = 0; i < dim; i++ )
  2113. -            sum += fabs(x->ve[i]);
  2114. -    else if ( scale->dim < dim )
  2115. -        error(E_SIZES,"_v_norm1");
  2116. -    else
  2117. -        for ( i = 0; i < dim; i++ )
  2118. -        {    s = scale->ve[i];
  2119. -            sum += ( s== 0.0 ) ? fabs(x->ve[i]) : fabs(x->ve[i]/s);
  2120. -        }
  2121. -
  2122. -    return sum;
  2123. -}
  2124. -
  2125. -/* square -- returns x^2 */
  2126. -double    square(x)
  2127. -double    x;
  2128. -{    return x*x;    }
  2129. -
  2130. -/* cube -- returns x^3 */
  2131. -double cube(x)
  2132. -double x;
  2133. -{  return x*x*x;   }
  2134. -
  2135. -/* _v_norm2 -- computes (scaled) 2-norm (Euclidean norm) of vectors */
  2136. -double    _v_norm2(x,scale)
  2137. -VEC    *x, *scale;
  2138. -{
  2139. -    int    i, dim;
  2140. -    Real    s, sum;
  2141. -
  2142. -    if ( x == (VEC *)NULL )
  2143. -        error(E_NULL,"_v_norm2");
  2144. -    dim = x->dim;
  2145. -
  2146. -    sum = 0.0;
  2147. -    if ( scale == (VEC *)NULL )
  2148. -        for ( i = 0; i < dim; i++ )
  2149. -            sum += square(x->ve[i]);
  2150. -    else if ( scale->dim < dim )
  2151. -        error(E_SIZES,"_v_norm2");
  2152. -    else
  2153. -        for ( i = 0; i < dim; i++ )
  2154. -        {    s = scale->ve[i];
  2155. -            sum += ( s== 0.0 ) ? square(x->ve[i]) :
  2156. -                            square(x->ve[i]/s);
  2157. -        }
  2158. -
  2159. -    return sqrt(sum);
  2160. -}
  2161. -
  2162. -#define    max(a,b)    ((a) > (b) ? (a) : (b))
  2163. -
  2164. -/* _v_norm_inf -- computes (scaled) infinity-norm (supremum norm) of vectors */
  2165. -double    _v_norm_inf(x,scale)
  2166. -VEC    *x, *scale;
  2167. -{
  2168. -    int    i, dim;
  2169. -    Real    s, maxval, tmp;
  2170. -
  2171. -    if ( x == (VEC *)NULL )
  2172. -        error(E_NULL,"_v_norm_inf");
  2173. -    dim = x->dim;
  2174. -
  2175. -    maxval = 0.0;
  2176. -    if ( scale == (VEC *)NULL )
  2177. -        for ( i = 0; i < dim; i++ )
  2178. -        {    tmp = fabs(x->ve[i]);
  2179. -            maxval = max(maxval,tmp);
  2180. -        }
  2181. -    else if ( scale->dim < dim )
  2182. -        error(E_SIZES,"_v_norm_inf");
  2183. -    else
  2184. -        for ( i = 0; i < dim; i++ )
  2185. -        {    s = scale->ve[i];
  2186. -            tmp = ( s== 0.0 ) ? fabs(x->ve[i]) : fabs(x->ve[i]/s);
  2187. -            maxval = max(maxval,tmp);
  2188. -        }
  2189. -
  2190. -    return maxval;
  2191. -}
  2192. -
  2193. -/* m_norm1 -- compute matrix 1-norm -- unscaled */
  2194. -double    m_norm1(A)
  2195. -MAT    *A;
  2196. -{
  2197. -    int    i, j, m, n;
  2198. -    Real    maxval, sum;
  2199. -
  2200. -    if ( A == (MAT *)NULL )
  2201. -        error(E_NULL,"m_norm1");
  2202. -
  2203. -    m = A->m;    n = A->n;
  2204. -    maxval = 0.0;
  2205. -
  2206. -    for ( j = 0; j < n; j++ )
  2207. -    {
  2208. -        sum = 0.0;
  2209. -        for ( i = 0; i < m; i ++ )
  2210. -            sum += fabs(A->me[i][j]);
  2211. -        maxval = max(maxval,sum);
  2212. -    }
  2213. -
  2214. -    return maxval;
  2215. -}
  2216. -
  2217. -/* m_norm_inf -- compute matrix infinity-norm -- unscaled */
  2218. -double    m_norm_inf(A)
  2219. -MAT    *A;
  2220. -{
  2221. -    int    i, j, m, n;
  2222. -    Real    maxval, sum;
  2223. -
  2224. -    if ( A == (MAT *)NULL )
  2225. -        error(E_NULL,"m_norm_inf");
  2226. -
  2227. -    m = A->m;    n = A->n;
  2228. -    maxval = 0.0;
  2229. -
  2230. -    for ( i = 0; i < m; i++ )
  2231. -    {
  2232. -        sum = 0.0;
  2233. -        for ( j = 0; j < n; j ++ )
  2234. -            sum += fabs(A->me[i][j]);
  2235. -        maxval = max(maxval,sum);
  2236. -    }
  2237. -
  2238. -    return maxval;
  2239. -}
  2240. -
  2241. -/* m_norm_frob -- compute matrix frobenius-norm -- unscaled */
  2242. -double    m_norm_frob(A)
  2243. -MAT    *A;
  2244. -{
  2245. -    int    i, j, m, n;
  2246. -    Real    sum;
  2247. -
  2248. -    if ( A == (MAT *)NULL )
  2249. -        error(E_NULL,"m_norm_frob");
  2250. -
  2251. -    m = A->m;    n = A->n;
  2252. -    sum = 0.0;
  2253. -
  2254. -    for ( i = 0; i < m; i++ )
  2255. -        for ( j = 0; j < n; j ++ )
  2256. -            sum += square(A->me[i][j]);
  2257. -
  2258. -    return sqrt(sum);
  2259. -}
  2260. -
  2261. //GO.SYSIN DD norm.c
  2262. echo hessen.c 1>&2
  2263. sed >hessen.c <<'//GO.SYSIN DD hessen.c' 's/^-//'
  2264. -
  2265. -/**************************************************************************
  2266. -**
  2267. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2268. -**
  2269. -**                 Meschach Library
  2270. -** 
  2271. -** This Meschach Library is provided "as is" without any express 
  2272. -** or implied warranty of any kind with respect to this software. 
  2273. -** In particular the authors shall not be liable for any direct, 
  2274. -** indirect, special, incidental or consequential damages arising 
  2275. -** in any way from use of the software.
  2276. -** 
  2277. -** Everyone is granted permission to copy, modify and redistribute this
  2278. -** Meschach Library, provided:
  2279. -**  1.  All copies contain this copyright notice.
  2280. -**  2.  All modified copies shall carry a notice stating who
  2281. -**      made the last modification and the date of such modification.
  2282. -**  3.  No charge is made for this software or works derived from it.  
  2283. -**      This clause shall not be construed as constraining other software
  2284. -**      distributed on the same medium as this software, nor is a
  2285. -**      distribution fee considered a charge.
  2286. -**
  2287. -***************************************************************************/
  2288. -
  2289. -
  2290. -
  2291. -/*
  2292. -        File containing routines for determining Hessenberg
  2293. -    factorisations.
  2294. -*/
  2295. -
  2296. -static    char    rcsid[] = "$Id: hessen.c,v 1.2 1994/01/13 05:36:24 des Exp $";
  2297. -
  2298. -#include    <stdio.h>
  2299. -#include    "matrix.h"
  2300. -#include        "matrix2.h"
  2301. -
  2302. -
  2303. -
  2304. -/* Hfactor -- compute Hessenberg factorisation in compact form.
  2305. -    -- factorisation performed in situ
  2306. -    -- for details of the compact form see QRfactor.c and matrix2.doc */
  2307. -MAT    *Hfactor(A, diag, beta)
  2308. -MAT    *A;
  2309. -VEC    *diag, *beta;
  2310. -{
  2311. -    static    VEC    *tmp1 = VNULL;
  2312. -    int    k, limit;
  2313. -
  2314. -    if ( ! A || ! diag || ! beta )
  2315. -        error(E_NULL,"Hfactor");
  2316. -    if ( diag->dim < A->m - 1 || beta->dim < A->m - 1 )
  2317. -        error(E_SIZES,"Hfactor");
  2318. -    if ( A->m != A->n )
  2319. -        error(E_SQUARE,"Hfactor");
  2320. -    limit = A->m - 1;
  2321. -
  2322. -    tmp1 = v_resize(tmp1,A->m);
  2323. -    MEM_STAT_REG(tmp1,TYPE_VEC);
  2324. -
  2325. -    for ( k = 0; k < limit; k++ )
  2326. -    {
  2327. -        get_col(A,(u_int)k,tmp1);
  2328. -        /* printf("the %d'th column = ");    v_output(tmp1); */
  2329. -        hhvec(tmp1,k+1,&beta->ve[k],tmp1,&A->me[k+1][k]);
  2330. -        /* diag->ve[k] = tmp1->ve[k+1]; */
  2331. -        v_set_val(diag,k,v_entry(tmp1,k+1));
  2332. -        /* printf("H/h vector = ");    v_output(tmp1); */
  2333. -        /* printf("from the %d'th entry\n",k+1); */
  2334. -        /* printf("beta = %g\n",beta->ve[k]); */
  2335. -
  2336. -        /* hhtrcols(A,k+1,k+1,tmp1,beta->ve[k]); */
  2337. -        /* hhtrrows(A,0  ,k+1,tmp1,beta->ve[k]); */
  2338. -        hhtrcols(A,k+1,k+1,tmp1,v_entry(beta,k));
  2339. -        hhtrrows(A,0  ,k+1,tmp1,v_entry(beta,k));
  2340. -        /* printf("A = ");        m_output(A); */
  2341. -    }
  2342. -
  2343. -    return (A);
  2344. -}
  2345. -
  2346. -/* makeHQ -- construct the Hessenberg orthogonalising matrix Q;
  2347. -    -- i.e. Hess M = Q.M.Q'    */
  2348. -MAT    *makeHQ(H, diag, beta, Qout)
  2349. -MAT    *H, *Qout;
  2350. -VEC    *diag, *beta;
  2351. -{
  2352. -    int    i, j, limit;
  2353. -    static    VEC    *tmp1 = VNULL, *tmp2 = VNULL;
  2354. -
  2355. -    if ( H==(MAT *)NULL || diag==(VEC *)NULL || beta==(VEC *)NULL )
  2356. -        error(E_NULL,"makeHQ");
  2357. -    limit = H->m - 1;
  2358. -    if ( diag->dim < limit || beta->dim < limit )
  2359. -        error(E_SIZES,"makeHQ");
  2360. -    if ( H->m != H->n )
  2361. -        error(E_SQUARE,"makeHQ");
  2362. -    Qout = m_resize(Qout,H->m,H->m);
  2363. -
  2364. -    tmp1 = v_resize(tmp1,H->m);
  2365. -    tmp2 = v_resize(tmp2,H->m);
  2366. -    MEM_STAT_REG(tmp1,TYPE_VEC);
  2367. -    MEM_STAT_REG(tmp2,TYPE_VEC);
  2368. -
  2369. -    for ( i = 0; i < H->m; i++ )
  2370. -    {
  2371. -        /* tmp1 = i'th basis vector */
  2372. -        for ( j = 0; j < H->m; j++ )
  2373. -            /* tmp1->ve[j] = 0.0; */
  2374. -            v_set_val(tmp1,j,0.0);
  2375. -        /* tmp1->ve[i] = 1.0; */
  2376. -        v_set_val(tmp1,i,1.0);
  2377. -
  2378. -        /* apply H/h transforms in reverse order */
  2379. -        for ( j = limit-1; j >= 0; j-- )
  2380. -        {
  2381. -            get_col(H,(u_int)j,tmp2);
  2382. -            /* tmp2->ve[j+1] = diag->ve[j]; */
  2383. -            v_set_val(tmp2,j+1,v_entry(diag,j));
  2384. -            hhtrvec(tmp2,beta->ve[j],j+1,tmp1,tmp1);
  2385. -        }
  2386. -
  2387. -        /* insert into Qout */
  2388. -        set_col(Qout,(u_int)i,tmp1);
  2389. -    }
  2390. -
  2391. -    return (Qout);
  2392. -}
  2393. -
  2394. -/* makeH -- construct actual Hessenberg matrix */
  2395. -MAT    *makeH(H,Hout)
  2396. -MAT    *H, *Hout;
  2397. -{
  2398. -    int    i, j, limit;
  2399. -
  2400. -    if ( H==(MAT *)NULL )
  2401. -        error(E_NULL,"makeH");
  2402. -    if ( H->m != H->n )
  2403. -        error(E_SQUARE,"makeH");
  2404. -    Hout = m_resize(Hout,H->m,H->m);
  2405. -    Hout = m_copy(H,Hout);
  2406. -
  2407. -    limit = H->m;
  2408. -    for ( i = 1; i < limit; i++ )
  2409. -        for ( j = 0; j < i-1; j++ )
  2410. -            /* Hout->me[i][j] = 0.0;*/
  2411. -            m_set_val(Hout,i,j,0.0);
  2412. -
  2413. -    return (Hout);
  2414. -}
  2415. -
  2416. //GO.SYSIN DD hessen.c
  2417. echo symmeig.c 1>&2
  2418. sed >symmeig.c <<'//GO.SYSIN DD symmeig.c' 's/^-//'
  2419. -
  2420. -/**************************************************************************
  2421. -**
  2422. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2423. -**
  2424. -**                 Meschach Library
  2425. -** 
  2426. -** This Meschach Library is provided "as is" without any express 
  2427. -** or implied warranty of any kind with respect to this software. 
  2428. -** In particular the authors shall not be liable for any direct, 
  2429. -** indirect, special, incidental or consequential damages arising 
  2430. -** in any way from use of the software.
  2431. -** 
  2432. -** Everyone is granted permission to copy, modify and redistribute this
  2433. -** Meschach Library, provided:
  2434. -**  1.  All copies contain this copyright notice.
  2435. -**  2.  All modified copies shall carry a notice stating who
  2436. -**      made the last modification and the date of such modification.
  2437. -**  3.  No charge is made for this software or works derived from it.  
  2438. -**      This clause shall not be construed as constraining other software
  2439. -**      distributed on the same medium as this software, nor is a
  2440. -**      distribution fee considered a charge.
  2441. -**
  2442. -***************************************************************************/
  2443. -
  2444. -
  2445. -/*
  2446. -    File containing routines for symmetric eigenvalue problems
  2447. -*/
  2448. -
  2449. -#include    <stdio.h>
  2450. -#include    <math.h>
  2451. -#include    "matrix.h"
  2452. -#include        "matrix2.h"
  2453. -
  2454. -
  2455. -static char rcsid[] = "$Id: symmeig.c,v 1.5 1994/02/16 03:23:39 des Exp $";
  2456. -
  2457. -
  2458. -
  2459. -#define    SQRT2    1.4142135623730949
  2460. -#define    sgn(x)    ( (x) >= 0 ? 1 : -1 )
  2461. -
  2462. -/* trieig -- finds eigenvalues of symmetric tridiagonal matrices
  2463. -    -- matrix represented by a pair of vectors a (diag entries)
  2464. -        and b (sub- & super-diag entries)
  2465. -    -- eigenvalues in a on return */
  2466. -VEC    *trieig(a,b,Q)
  2467. -VEC    *a, *b;
  2468. -MAT    *Q;
  2469. -{
  2470. -    int    i, i_min, i_max, n, split;
  2471. -    Real    *a_ve, *b_ve;
  2472. -    Real    b_sqr, bk, ak1, bk1, ak2, bk2, z;
  2473. -    Real    c, c2, cs, s, s2, d, mu;
  2474. -
  2475. -    if ( ! a || ! b )
  2476. -        error(E_NULL,"trieig");
  2477. -    if ( a->dim != b->dim + 1 || ( Q && Q->m != a->dim ) )
  2478. -        error(E_SIZES,"trieig");
  2479. -    if ( Q && Q->m != Q->n )
  2480. -        error(E_SQUARE,"trieig");
  2481. -
  2482. -    n = a->dim;
  2483. -    a_ve = a->ve;        b_ve = b->ve;
  2484. -
  2485. -    i_min = 0;
  2486. -    while ( i_min < n )        /* outer while loop */
  2487. -    {
  2488. -        /* find i_max to suit;
  2489. -            submatrix i_min..i_max should be irreducible */
  2490. -        i_max = n-1;
  2491. -        for ( i = i_min; i < n-1; i++ )
  2492. -            if ( b_ve[i] == 0.0 )
  2493. -            {    i_max = i;    break;    }
  2494. -        if ( i_max <= i_min )
  2495. -        {
  2496. -            /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  2497. -            i_min = i_max + 1;
  2498. -            continue;    /* outer while loop */
  2499. -        }
  2500. -
  2501. -        /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  2502. -
  2503. -        /* repeatedly perform QR method until matrix splits */
  2504. -        split = FALSE;
  2505. -        while ( ! split )        /* inner while loop */
  2506. -        {
  2507. -
  2508. -            /* find Wilkinson shift */
  2509. -            d = (a_ve[i_max-1] - a_ve[i_max])/2;
  2510. -            b_sqr = b_ve[i_max-1]*b_ve[i_max-1];
  2511. -            mu = a_ve[i_max] - b_sqr/(d + sgn(d)*sqrt(d*d+b_sqr));
  2512. -            /* printf("# Wilkinson shift = %g\n",mu); */
  2513. -
  2514. -            /* initial Givens' rotation */
  2515. -            givens(a_ve[i_min]-mu,b_ve[i_min],&c,&s);
  2516. -            s = -s;
  2517. -            /* printf("# c = %g, s = %g\n",c,s); */
  2518. -            if ( fabs(c) < SQRT2 )
  2519. -            {    c2 = c*c;    s2 = 1-c2;    }
  2520. -            else
  2521. -            {    s2 = s*s;    c2 = 1-s2;    }
  2522. -            cs = c*s;
  2523. -            ak1 = c2*a_ve[i_min]+s2*a_ve[i_min+1]-2*cs*b_ve[i_min];
  2524. -            bk1 = cs*(a_ve[i_min]-a_ve[i_min+1]) +
  2525. -                        (c2-s2)*b_ve[i_min];
  2526. -            ak2 = s2*a_ve[i_min]+c2*a_ve[i_min+1]+2*cs*b_ve[i_min];
  2527. -            bk2 = ( i_min < i_max-1 ) ? c*b_ve[i_min+1] : 0.0;
  2528. -            z  = ( i_min < i_max-1 ) ? -s*b_ve[i_min+1] : 0.0;
  2529. -            a_ve[i_min] = ak1;
  2530. -            a_ve[i_min+1] = ak2;
  2531. -            b_ve[i_min] = bk1;
  2532. -            if ( i_min < i_max-1 )
  2533. -            b_ve[i_min+1] = bk2;
  2534. -            if ( Q )
  2535. -            rot_cols(Q,i_min,i_min+1,c,-s,Q);
  2536. -            /* printf("# z = %g\n",z); */
  2537. -            /* printf("# a [temp1] =\n");    v_output(a); */
  2538. -            /* printf("# b [temp1] =\n");    v_output(b); */
  2539. -
  2540. -            for ( i = i_min+1; i < i_max; i++ )
  2541. -            {
  2542. -            /* get Givens' rotation for sub-block -- k == i-1 */
  2543. -            givens(b_ve[i-1],z,&c,&s);
  2544. -            s = -s;
  2545. -            /* printf("# c = %g, s = %g\n",c,s); */
  2546. -
  2547. -            /* perform Givens' rotation on sub-block */
  2548. -                if ( fabs(c) < SQRT2 )
  2549. -                {    c2 = c*c;    s2 = 1-c2;    }
  2550. -                else
  2551. -                {    s2 = s*s;    c2 = 1-s2;    }
  2552. -                cs = c*s;
  2553. -            bk  = c*b_ve[i-1] - s*z;
  2554. -            ak1 = c2*a_ve[i]+s2*a_ve[i+1]-2*cs*b_ve[i];
  2555. -            bk1 = cs*(a_ve[i]-a_ve[i+1]) +
  2556. -                        (c2-s2)*b_ve[i];
  2557. -            ak2 = s2*a_ve[i]+c2*a_ve[i+1]+2*cs*b_ve[i];
  2558. -            bk2 = ( i+1 < i_max ) ? c*b_ve[i+1] : 0.0;
  2559. -            z  = ( i+1 < i_max ) ? -s*b_ve[i+1] : 0.0;
  2560. -            a_ve[i] = ak1;    a_ve[i+1] = ak2;
  2561. -            b_ve[i] = bk1;
  2562. -            if ( i < i_max-1 )
  2563. -                b_ve[i+1] = bk2;
  2564. -            if ( i > i_min )
  2565. -                b_ve[i-1] = bk;
  2566. -            if ( Q )
  2567. -                rot_cols(Q,i,i+1,c,-s,Q);
  2568. -                /* printf("# a [temp2] =\n");    v_output(a); */
  2569. -                /* printf("# b [temp2] =\n");    v_output(b); */
  2570. -            }
  2571. -
  2572. -            /* test to see if matrix should be split */
  2573. -            for ( i = i_min; i < i_max; i++ )
  2574. -            if ( fabs(b_ve[i]) < MACHEPS*
  2575. -                    (fabs(a_ve[i])+fabs(a_ve[i+1])) )
  2576. -            {   b_ve[i] = 0.0;    split = TRUE;    }
  2577. -
  2578. -            /* printf("# a =\n");    v_output(a); */
  2579. -            /* printf("# b =\n");    v_output(b); */
  2580. -        }
  2581. -    }
  2582. -
  2583. -    return a;
  2584. -}
  2585. -
  2586. -/* symmeig -- computes eigenvalues of a dense symmetric matrix
  2587. -    -- A **must** be symmetric on entry
  2588. -    -- eigenvalues stored in out
  2589. -    -- Q contains orthogonal matrix of eigenvectors
  2590. -    -- returns vector of eigenvalues */
  2591. -VEC    *symmeig(A,Q,out)
  2592. -MAT    *A, *Q;
  2593. -VEC    *out;
  2594. -{
  2595. -    int    i;
  2596. -    static MAT    *tmp = MNULL;
  2597. -    static VEC    *b   = VNULL, *diag = VNULL, *beta = VNULL;
  2598. -
  2599. -    if ( ! A )
  2600. -        error(E_NULL,"symmeig");
  2601. -    if ( A->m != A->n )
  2602. -        error(E_SQUARE,"symmeig");
  2603. -    if ( ! out || out->dim != A->m )
  2604. -        out = v_resize(out,A->m);
  2605. -
  2606. -    tmp  = m_copy(A,tmp);
  2607. -    b    = v_resize(b,A->m - 1);
  2608. -    diag = v_resize(diag,(u_int)A->m);
  2609. -    beta = v_resize(beta,(u_int)A->m);
  2610. -    MEM_STAT_REG(tmp,TYPE_MAT);
  2611. -    MEM_STAT_REG(b,TYPE_VEC);
  2612. -    MEM_STAT_REG(diag,TYPE_VEC);
  2613. -    MEM_STAT_REG(beta,TYPE_VEC);
  2614. -
  2615. -    Hfactor(tmp,diag,beta);
  2616. -    if ( Q )
  2617. -        makeHQ(tmp,diag,beta,Q);
  2618. -
  2619. -    for ( i = 0; i < A->m - 1; i++ )
  2620. -    {
  2621. -        out->ve[i] = tmp->me[i][i];
  2622. -        b->ve[i] = tmp->me[i][i+1];
  2623. -    }
  2624. -    out->ve[i] = tmp->me[i][i];
  2625. -    trieig(out,b,Q);
  2626. -
  2627. -    return out;
  2628. -}
  2629. -
  2630. //GO.SYSIN DD symmeig.c
  2631. echo schur.c 1>&2
  2632. sed >schur.c <<'//GO.SYSIN DD schur.c' 's/^-//'
  2633. -
  2634. -/**************************************************************************
  2635. -**
  2636. -** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
  2637. -**
  2638. -**                 Meschach Library
  2639. -** 
  2640. -** This Meschach Library is provided "as is" without any express 
  2641. -** or implied warranty of any kind with respect to this software. 
  2642. -** In particular the authors shall not be liable for any direct, 
  2643. -** indirect, special, incidental or consequential damages arising 
  2644. -** in any way from use of the software.
  2645. -** 
  2646. -** Everyone is granted permission to copy, modify and redistribute this
  2647. -** Meschach Library, provided:
  2648. -**  1.  All copies contain this copyright notice.
  2649. -**  2.  All modified copies shall carry a notice stating who
  2650. -**      made the last modification and the date of such modification.
  2651. -**  3.  No charge is made for this software or works derived from it.  
  2652. -**      This clause shall not be construed as constraining other software
  2653. -**      distributed on the same medium as this software, nor is a
  2654. -**      distribution fee considered a charge.
  2655. -**
  2656. -***************************************************************************/
  2657. -
  2658. -
  2659. -/*    
  2660. -    File containing routines for computing the Schur decomposition
  2661. -    of a real non-symmetric matrix
  2662. -    See also: hessen.c
  2663. -*/
  2664. -
  2665. -#include    <stdio.h>
  2666. -#include    <math.h>
  2667. -#include    "matrix.h"
  2668. -#include        "matrix2.h"
  2669. -
  2670. -
  2671. -static char rcsid[] = "$Id: schur.c,v 1.7 1994/03/17 05:36:53 des Exp $";
  2672. -
  2673. -
  2674. -
  2675. -#ifndef ANSI_C
  2676. -static    void    hhldr3(x,y,z,nu1,beta,newval)
  2677. -double    x, y, z;
  2678. -Real    *nu1, *beta, *newval;
  2679. -#else
  2680. -static    void    hhldr3(double x, double y, double z,
  2681. -               Real *nu1, Real *beta, Real *newval)
  2682. -#endif
  2683. -{
  2684. -    Real    alpha;
  2685. -
  2686. -    if ( x >= 0.0 )
  2687. -        alpha = sqrt(x*x+y*y+z*z);
  2688. -    else
  2689. -        alpha = -sqrt(x*x+y*y+z*z);
  2690. -    *nu1 = x + alpha;
  2691. -    *beta = 1.0/(alpha*(*nu1));
  2692. -    *newval = alpha;
  2693. -}
  2694. -
  2695. -#ifndef ANSI_C
  2696. -static    void    hhldr3cols(A,k,j0,beta,nu1,nu2,nu3)
  2697. -MAT    *A;
  2698. -int    k, j0;
  2699. -double    beta, nu1, nu2, nu3;
  2700. -#else
  2701. -static    void    hhldr3cols(MAT *A, int k, int j0, double beta,
  2702. -               double nu1, double nu2, double nu3)
  2703. -#endif
  2704. -{
  2705. -    Real    **A_me, ip, prod;
  2706. -    int    j, n;
  2707. -
  2708. -    if ( k < 0 || k+3 > A->m || j0 < 0 )
  2709. -        error(E_BOUNDS,"hhldr3cols");
  2710. -    A_me = A->me;        n = A->n;
  2711. -
  2712. -    /* printf("hhldr3cols:(l.%d) j0 = %d, k = %d, A at 0x%lx, m = %d, n = %d\n",
  2713. -           __LINE__, j0, k, (long)A, A->m, A->n); */
  2714. -    /* printf("hhldr3cols: A (dumped) =\n");    m_dump(stdout,A); */
  2715. -
  2716. -    for ( j = j0; j < n; j++ )
  2717. -    {
  2718. -        /*****        
  2719. -        ip = nu1*A_me[k][j] + nu2*A_me[k+1][j] + nu3*A_me[k+2][j];
  2720. -        prod = ip*beta;
  2721. -        A_me[k][j]   -= prod*nu1;
  2722. -        A_me[k+1][j] -= prod*nu2;
  2723. -        A_me[k+2][j] -= prod*nu3;
  2724. -        *****/
  2725. -        /* printf("hhldr3cols: j = %d\n", j); */
  2726. -
  2727. -        ip = nu1*m_entry(A,k,j)+nu2*m_entry(A,k+1,j)+nu3*m_entry(A,k+2,j);
  2728. -        prod = ip*beta;
  2729. -        /*****
  2730. -        m_set_val(A,k  ,j,m_entry(A,k  ,j) - prod*nu1);
  2731. -        m_set_val(A,k+1,j,m_entry(A,k+1,j) - prod*nu2);
  2732. -        m_set_val(A,k+2,j,m_entry(A,k+2,j) - prod*nu3);
  2733. -        *****/
  2734. -        m_add_val(A,k  ,j,-prod*nu1);
  2735. -        m_add_val(A,k+1,j,-prod*nu2);
  2736. -        m_add_val(A,k+2,j,-prod*nu3);
  2737. -
  2738. -    }
  2739. -    /* printf("hhldr3cols:(l.%d) j0 = %d, k = %d, m = %d, n = %d\n",
  2740. -           __LINE__, j0, k, A->m, A->n); */
  2741. -    /* putc('\n',stdout); */
  2742. -}
  2743. -
  2744. -#ifndef ANSI_C
  2745. -static    void    hhldr3rows(A,k,i0,beta,nu1,nu2,nu3)
  2746. -MAT    *A;
  2747. -int    k, i0;
  2748. -double    beta, nu1, nu2, nu3;
  2749. -#else
  2750. -static    void    hhldr3rows(MAT *A, int k, int i0, double beta,
  2751. -               double nu1, double nu2, double nu3)
  2752. -#endif
  2753. -{
  2754. -    Real    **A_me, ip, prod;
  2755. -    int    i, m;
  2756. -
  2757. -    /* printf("hhldr3rows:(l.%d) A at 0x%lx\n", __LINE__, (long)A); */
  2758. -    /* printf("hhldr3rows: k = %d\n", k); */
  2759. -    if ( k < 0 || k+3 > A->n )
  2760. -        error(E_BOUNDS,"hhldr3rows");
  2761. -    A_me = A->me;        m = A->m;
  2762. -    i0 = min(i0,m-1);
  2763. -
  2764. -    for ( i = 0; i <= i0; i++ )
  2765. -    {
  2766. -        /****
  2767. -        ip = nu1*A_me[i][k] + nu2*A_me[i][k+1] + nu3*A_me[i][k+2];
  2768. -        prod = ip*beta;
  2769. -        A_me[i][k]   -= prod*nu1;
  2770. -        A_me[i][k+1] -= prod*nu2;
  2771. -        A_me[i][k+2] -= prod*nu3;
  2772. -        ****/
  2773. -
  2774. -        ip = nu1*m_entry(A,i,k)+nu2*m_entry(A,i,k+1)+nu3*m_entry(A,i,k+2);
  2775. -        prod = ip*beta;
  2776. -        m_add_val(A,i,k  , - prod*nu1);
  2777. -        m_add_val(A,i,k+1, - prod*nu2);
  2778. -        m_add_val(A,i,k+2, - prod*nu3);
  2779. -
  2780. -    }
  2781. -}
  2782. -
  2783. -/* schur -- computes the Schur decomposition of the matrix A in situ
  2784. -    -- optionally, gives Q matrix such that Q^T.A.Q is upper triangular
  2785. -    -- returns upper triangular Schur matrix */
  2786. -MAT    *schur(A,Q)
  2787. -MAT    *A, *Q;
  2788. -{
  2789. -    int        i, j, iter, k, k_min, k_max, k_tmp, n, split;
  2790. -    Real    beta2, c, discrim, dummy, nu1, s, t, tmp, x, y, z;
  2791. -    Real    **A_me;
  2792. -    Real    sqrt_macheps;
  2793. -    static    VEC    *diag=VNULL, *beta=VNULL;
  2794. -    
  2795. -    if ( ! A )
  2796. -    error(E_NULL,"schur");
  2797. -    if ( A->m != A->n || ( Q && Q->m != Q->n ) )
  2798. -    error(E_SQUARE,"schur");
  2799. -    if ( Q != MNULL && Q->m != A->m )
  2800. -    error(E_SIZES,"schur");
  2801. -    n = A->n;
  2802. -    diag = v_resize(diag,A->n);
  2803. -    beta = v_resize(beta,A->n);
  2804. -    MEM_STAT_REG(diag,TYPE_VEC);
  2805. -    MEM_STAT_REG(beta,TYPE_VEC);
  2806. -    /* compute Hessenberg form */
  2807. -    Hfactor(A,diag,beta);
  2808. -    
  2809. -    /* save Q if necessary */
  2810. -    if ( Q )
  2811. -    Q = makeHQ(A,diag,beta,Q);
  2812. -    makeH(A,A);
  2813. -
  2814. -    sqrt_macheps = sqrt(MACHEPS);
  2815. -
  2816. -    k_min = 0;    A_me = A->me;
  2817. -
  2818. -    while ( k_min < n )
  2819. -    {
  2820. -    Real    a00, a01, a10, a11;
  2821. -    double    scale, t, numer, denom;
  2822. -
  2823. -    /* find k_max to suit:
  2824. -       submatrix k_min..k_max should be irreducible */
  2825. -    k_max = n-1;
  2826. -    for ( k = k_min; k < k_max; k++ )
  2827. -        /* if ( A_me[k+1][k] == 0.0 ) */
  2828. -        if ( m_entry(A,k+1,k) == 0.0 )
  2829. -        {    k_max = k;    break;    }
  2830. -
  2831. -    if ( k_max <= k_min )
  2832. -    {
  2833. -        k_min = k_max + 1;
  2834. -        continue;        /* outer loop */
  2835. -    }
  2836. -
  2837. -    /* check to see if we have a 2 x 2 block
  2838. -       with complex eigenvalues */
  2839. -    if ( k_max == k_min + 1 )
  2840. -    {
  2841. -        /* tmp = A_me[k_min][k_min] - A_me[k_max][k_max]; */
  2842. -        a00 = m_entry(A,k_min,k_min);
  2843. -        a01 = m_entry(A,k_min,k_max);
  2844. -        a10 = m_entry(A,k_max,k_min);
  2845. -        a11 = m_entry(A,k_max,k_max);
  2846. -        tmp = a00 - a11;
  2847. -        /* discrim = tmp*tmp +
  2848. -        4*A_me[k_min][k_max]*A_me[k_max][k_min]; */
  2849. -        discrim = tmp*tmp +
  2850. -        4*a01*a10;
  2851. -        if ( discrim < 0.0 )
  2852. -        {    /* yes -- e-vals are complex
  2853. -           -- put 2 x 2 block in form [a b; c a];
  2854. -           then eigenvalues have real part a & imag part sqrt(|bc|) */
  2855. -        numer = - tmp;
  2856. -        denom = ( a01+a10 >= 0.0 ) ?
  2857. -            (a01+a10) + sqrt((a01+a10)*(a01+a10)+tmp*tmp) :
  2858. -            (a01+a10) - sqrt((a01+a10)*(a01+a10)+tmp*tmp);
  2859. -        if ( denom != 0.0 )
  2860. -        {   /* t = s/c = numer/denom */
  2861. -            t = numer/denom;
  2862. -            scale = c = 1.0/sqrt(1+t*t);
  2863. -            s = c*t;
  2864. -        }
  2865. -        else
  2866. -        {
  2867. -            c = 1.0;
  2868. -            s = 0.0;
  2869. -        }
  2870. -        rot_cols(A,k_min,k_max,c,s,A);
  2871. -        rot_rows(A,k_min,k_max,c,s,A);
  2872. -        if ( Q != MNULL )
  2873. -            rot_cols(Q,k_min,k_max,c,s,Q);
  2874. -        k_min = k_max + 1;
  2875. -        continue;
  2876. -        }
  2877. -        else /* discrim >= 0; i.e. block has two real eigenvalues */
  2878. -        {    /* no -- e-vals are not complex;
  2879. -           split 2 x 2 block and continue */
  2880. -        /* s/c = numer/denom */
  2881. -        numer = ( tmp >= 0.0 ) ?
  2882. -            - tmp - sqrt(discrim) : - tmp + sqrt(discrim);
  2883. -        denom = 2*a01;
  2884. -        if ( fabs(numer) < fabs(denom) )
  2885. -        {   /* t = s/c = numer/denom */
  2886. -            t = numer/denom;
  2887. -            scale = c = 1.0/sqrt(1+t*t);
  2888. -            s = c*t;
  2889. -        }
  2890. -        else if ( numer != 0.0 )
  2891. -        {   /* t = c/s = denom/numer */
  2892. -            t = denom/numer;
  2893. -            scale = 1.0/sqrt(1+t*t);
  2894. -            c = fabs(t)*scale;
  2895. -            s = ( t >= 0.0 ) ? scale : -scale;
  2896. -        }
  2897. -        else /* numer == denom == 0 */
  2898. -        {
  2899. -            c = 0.0;
  2900. -            s = 1.0;
  2901. -        }
  2902. -        rot_cols(A,k_min,k_max,c,s,A);
  2903. -        rot_rows(A,k_min,k_max,c,s,A);
  2904. -        /* A->me[k_max][k_min] = 0.0; */
  2905. -        if ( Q != MNULL )
  2906. -            rot_cols(Q,k_min,k_max,c,s,Q);
  2907. -        k_min = k_max + 1;    /* go to next block */
  2908. -        continue;
  2909. -        }
  2910. -    }
  2911. -
  2912. -    /* now have r x r block with r >= 2:
  2913. -       apply Francis QR step until block splits */
  2914. -    split = FALSE;        iter = 0;
  2915. -    while ( ! split )
  2916. -    {
  2917. -        iter++;
  2918. -        
  2919. -        /* set up Wilkinson/Francis complex shift */
  2920. -        k_tmp = k_max - 1;
  2921. -
  2922. -        a00 = m_entry(A,k_tmp,k_tmp);
  2923. -        a01 = m_entry(A,k_tmp,k_max);
  2924. -        a10 = m_entry(A,k_max,k_tmp);
  2925. -        a11 = m_entry(A,k_max,k_max);
  2926. -
  2927. -        /* treat degenerate cases differently
  2928. -           -- if there are still no splits after five iterations
  2929. -              and the bottom 2 x 2 looks degenerate, force it to
  2930. -          split */
  2931. -        if ( iter >= 5 &&
  2932. -         fabs(a00-a11) < sqrt_macheps*(fabs(a00)+fabs(a11)) &&
  2933. -         (fabs(a01) < sqrt_macheps*(fabs(a00)+fabs(a11)) ||
  2934. -          fabs(a10) < sqrt_macheps*(fabs(a00)+fabs(a11))) )
  2935. -        {
  2936. -          if ( fabs(a01) < sqrt_macheps*(fabs(a00)+fabs(a11)) )
  2937. -        m_set_val(A,k_tmp,k_max,0.0);
  2938. -          if ( fabs(a10) < sqrt_macheps*(fabs(a00)+fabs(a11)) )
  2939. -        {
  2940. -          m_set_val(A,k_max,k_tmp,0.0);
  2941. -          split = TRUE;
  2942. -          continue;
  2943. -        }
  2944. -        }
  2945. -
  2946. -        s = a00 + a11;
  2947. -        t = a00*a11 - a01*a10;
  2948. -
  2949. -        /* break loop if a 2 x 2 complex block */
  2950. -        if ( k_max == k_min + 1 && s*s < 4.0*t )
  2951. -        {
  2952. -        split = TRUE;
  2953. -        continue;
  2954. -        }
  2955. -
  2956. -        /* perturb shift if convergence is slow */
  2957. -        if ( (iter % 10) == 0 )
  2958. -        {    s += iter*0.02;        t += iter*0.02;
  2959. -        }
  2960. -
  2961. -        /* set up Householder transformations */
  2962. -        k_tmp = k_min + 1;
  2963. -        /********************
  2964. -        x = A_me[k_min][k_min]*A_me[k_min][k_min] +
  2965. -        A_me[k_min][k_tmp]*A_me[k_tmp][k_min] -
  2966. -            s*A_me[k_min][k_min] + t;
  2967. -        y = A_me[k_tmp][k_min]*
  2968. -        (A_me[k_min][k_min]+A_me[k_tmp][k_tmp]-s);
  2969. -        if ( k_min + 2 <= k_max )
  2970. -        z = A_me[k_tmp][k_min]*A_me[k_min+2][k_tmp];
  2971. -        else
  2972. -        z = 0.0;
  2973. -        ********************/
  2974. -
  2975. -        a00 = m_entry(A,k_min,k_min);
  2976. -        a01 = m_entry(A,k_min,k_tmp);
  2977. -        a10 = m_entry(A,k_tmp,k_min);
  2978. -        a11 = m_entry(A,k_tmp,k_tmp);
  2979. -
  2980. -        /********************
  2981. -        a00 = A->me[k_min][k_min];
  2982. -        a01 = A->me[k_min][k_tmp];
  2983. -        a10 = A->me[k_tmp][k_min];
  2984. -        a11 = A->me[k_tmp][k_tmp];
  2985. -        ********************/
  2986. -        x = a00*a00 + a01*a10 - s*a00 + t;
  2987. -        y = a10*(a00+a11-s);
  2988. -        if ( k_min + 2 <= k_max )
  2989. -        z = a10* /* m_entry(A,k_min+2,k_tmp) */ A->me[k_min+2][k_tmp];
  2990. -        else
  2991. -        z = 0.0;
  2992. -
  2993. -        for ( k = k_min; k <= k_max-1; k++ )
  2994. -        {
  2995. -        if ( k < k_max - 1 )
  2996. -        {
  2997. -            hhldr3(x,y,z,&nu1,&beta2,&dummy);
  2998. -            tracecatch(hhldr3cols(A,k,max(k-1,0),  beta2,nu1,y,z),"schur");
  2999. -            tracecatch(hhldr3rows(A,k,min(n-1,k+3),beta2,nu1,y,z),"schur");
  3000. -            if ( Q != MNULL )
  3001. -            hhldr3rows(Q,k,n-1,beta2,nu1,y,z);
  3002. -        }
  3003. -        else
  3004. -        {
  3005. -            givens(x,y,&c,&s);
  3006. -            rot_cols(A,k,k+1,c,s,A);
  3007. -            rot_rows(A,k,k+1,c,s,A);
  3008. -            if ( Q )
  3009. -            rot_cols(Q,k,k+1,c,s,Q);
  3010. -        }
  3011. -        /* if ( k >= 2 )
  3012. -            m_set_val(A,k,k-2,0.0); */
  3013. -        /* x = A_me[k+1][k]; */
  3014. -        x = m_entry(A,k+1,k);
  3015. -        if ( k <= k_max - 2 )
  3016. -            /* y = A_me[k+2][k];*/
  3017. -            y = m_entry(A,k+2,k);
  3018. -        else
  3019. -            y = 0.0;
  3020. -        if ( k <= k_max - 3 )
  3021. -            /* z = A_me[k+3][k]; */
  3022. -            z = m_entry(A,k+3,k);
  3023. -        else
  3024. -            z = 0.0;
  3025. -        }
  3026. -        /* if ( k_min > 0 )
  3027. -        m_set_val(A,k_min,k_min-1,0.0);
  3028. -        if ( k_max < n - 1 )
  3029. -        m_set_val(A,k_max+1,k_max,0.0); */
  3030. -        for ( k = k_min; k <= k_max-2; k++ )
  3031. -        {
  3032. -        /* zero appropriate sub-diagonals */
  3033. -        m_set_val(A,k+2,k,0.0);
  3034. -        if ( k < k_max-2 )
  3035. -            m_set_val(A,k+3,k,0.0);
  3036. -        }
  3037. -
  3038. -        /* test to see if matrix should split */
  3039. -        for ( k = k_min; k < k_max; k++ )
  3040. -        if ( fabs(A_me[k+1][k]) < MACHEPS*
  3041. -            (fabs(A_me[k][k])+fabs(A_me[k+1][k+1])) )
  3042. -        {    A_me[k+1][k] = 0.0;    split = TRUE;    }
  3043. -    }
  3044. -    }
  3045. -    
  3046. -    /* polish up A by zeroing strictly lower triangular elements
  3047. -       and small sub-diagonal elements */
  3048. -    for ( i = 0; i < A->m; i++ )
  3049. -    for ( j = 0; j < i-1; j++ )
  3050. -        A_me[i][j] = 0.0;
  3051. -    for ( i = 0; i < A->m - 1; i++ )
  3052. -    if ( fabs(A_me[i+1][i]) < MACHEPS*
  3053. -        (fabs(A_me[i][i])+fabs(A_me[i+1][i+1])) )
  3054. -        A_me[i+1][i] = 0.0;
  3055. -
  3056. -    return A;
  3057. -}
  3058. -
  3059. -/* schur_vals -- compute real & imaginary parts of eigenvalues
  3060. -    -- assumes T contains a block upper triangular matrix
  3061. -        as produced by schur()
  3062. -    -- real parts stored in real_pt, imaginary parts in imag_pt */
  3063. -void    schur_evals(T,real_pt,imag_pt)
  3064. -MAT    *T;
  3065. -VEC    *real_pt, *imag_pt;
  3066. -{
  3067. -    int    i, n;
  3068. -    Real    discrim, **T_me;
  3069. -    Real    diff, sum, tmp;
  3070. -
  3071. -    if ( ! T || ! real_pt || ! imag_pt )
  3072. -        error(E_NULL,"schur_evals");
  3073. -    if ( T->m != T->n )
  3074. -        error(E_SQUARE,"schur_evals");
  3075. -    n = T->n;    T_me = T->me;
  3076. -    real_pt = v_resize(real_pt,(u_int)n);
  3077. -    imag_pt = v_resize(imag_pt,(u_int)n);
  3078. -
  3079. -    i = 0;
  3080. -    while ( i < n )
  3081. -    {
  3082. -        if ( i < n-1 && T_me[i+1][i] != 0.0 )
  3083. -        {   /* should be a complex eigenvalue */
  3084. -            sum  = 0.5*(T_me[i][i]+T_me[i+1][i+1]);
  3085. -            diff = 0.5*(T_me[i][i]-T_me[i+1][i+1]);
  3086. -            discrim = diff*diff + T_me[i][i+1]*T_me[i+1][i];
  3087. -            if ( discrim < 0.0 )
  3088. -            {    /* yes -- complex e-vals */
  3089. -            real_pt->ve[i] = real_pt->ve[i+1] = sum;
  3090. -            imag_pt->ve[i] = sqrt(-discrim);
  3091. -            imag_pt->ve[i+1] = - imag_pt->ve[i];
  3092. -            }
  3093. -            else
  3094. -            {    /* no -- actually both real */
  3095. -            tmp = sqrt(discrim);
  3096. -            real_pt->ve[i]   = sum + tmp;
  3097. -            real_pt->ve[i+1] = sum - tmp;
  3098. -            imag_pt->ve[i]   = imag_pt->ve[i+1] = 0.0;
  3099. -            }
  3100. -            i += 2;
  3101. -        }
  3102. -        else
  3103. -        {   /* real eigenvalue */
  3104. -            real_pt->ve[i] = T_me[i][i];
  3105. -            imag_pt->ve[i] = 0.0;
  3106. -            i++;
  3107. -        }
  3108. -    }
  3109. -}
  3110. -
  3111. -/* schur_vecs -- returns eigenvectors computed from the real Schur
  3112. -        decomposition of a matrix
  3113. -    -- T is the block upper triangular Schur matrix
  3114. -    -- Q is the orthognal matrix where A = Q.T.Q^T
  3115. -    -- if Q is null, the eigenvectors of T are returned
  3116. -    -- X_re is the real part of the matrix of eigenvectors,
  3117. -        and X_im is the imaginary part of the matrix.
  3118. -    -- X_re is returned */
  3119. -MAT    *schur_vecs(T,Q,X_re,X_im)
  3120. -MAT    *T, *Q, *X_re, *X_im;
  3121. -{
  3122. -    int    i, j, limit;
  3123. -    Real    t11_re, t11_im, t12, t21, t22_re, t22_im;
  3124. -    Real    l_re, l_im, det_re, det_im, invdet_re, invdet_im,
  3125. -        val1_re, val1_im, val2_re, val2_im,
  3126. -        tmp_val1_re, tmp_val1_im, tmp_val2_re, tmp_val2_im, **T_me;
  3127. -    Real    sum, diff, discrim, magdet, norm, scale;
  3128. -    static VEC    *tmp1_re=VNULL, *tmp1_im=VNULL,
  3129. -            *tmp2_re=VNULL, *tmp2_im=VNULL;
  3130. -
  3131. -    if ( ! T || ! X_re )
  3132. -        error(E_NULL,"schur_vecs");
  3133. -    if ( T->m != T->n || X_re->m != X_re->n ||
  3134. -        ( Q != MNULL && Q->m != Q->n ) ||
  3135. -        ( X_im != MNULL && X_im->m != X_im->n ) )
  3136. -        error(E_SQUARE,"schur_vecs");
  3137. -    if ( T->m != X_re->m ||
  3138. -        ( Q != MNULL && T->m != Q->m ) ||
  3139. -        ( X_im != MNULL && T->m != X_im->m ) )
  3140. -        error(E_SIZES,"schur_vecs");
  3141. -
  3142. -    tmp1_re = v_resize(tmp1_re,T->m);
  3143. -    tmp1_im = v_resize(tmp1_im,T->m);
  3144. -    tmp2_re = v_resize(tmp2_re,T->m);
  3145. -    tmp2_im = v_resize(tmp2_im,T->m);
  3146. -    MEM_STAT_REG(tmp1_re,TYPE_VEC);
  3147. -    MEM_STAT_REG(tmp1_im,TYPE_VEC);
  3148. -    MEM_STAT_REG(tmp2_re,TYPE_VEC);
  3149. -    MEM_STAT_REG(tmp2_im,TYPE_VEC);
  3150. -
  3151. -    T_me = T->me;
  3152. -    i = 0;
  3153. -    while ( i < T->m )
  3154. -    {
  3155. -        if ( i+1 < T->m && T->me[i+1][i] != 0.0 )
  3156. -        {    /* complex eigenvalue */
  3157. -        sum  = 0.5*(T_me[i][i]+T_me[i+1][i+1]);
  3158. -        diff = 0.5*(T_me[i][i]-T_me[i+1][i+1]);
  3159. -        discrim = diff*diff + T_me[i][i+1]*T_me[i+1][i];
  3160. -        l_re = l_im = 0.0;
  3161. -        if ( discrim < 0.0 )
  3162. -        {    /* yes -- complex e-vals */
  3163. -            l_re = sum;
  3164. -            l_im = sqrt(-discrim);
  3165. -        }
  3166. -        else /* not correct Real Schur form */
  3167. -            error(E_RANGE,"schur_vecs");
  3168. -        }
  3169. -        else
  3170. -        {
  3171. -        l_re = T_me[i][i];
  3172. -        l_im = 0.0;
  3173. -        }
  3174. -
  3175. -        v_zero(tmp1_im);
  3176. -        v_rand(tmp1_re);
  3177. -        sv_mlt(MACHEPS,tmp1_re,tmp1_re);
  3178. -
  3179. -        /* solve (T-l.I)x = tmp1 */
  3180. -        limit = ( l_im != 0.0 ) ? i+1 : i;
  3181. -        /* printf("limit = %d\n",limit); */
  3182. -        for ( j = limit+1; j < T->m; j++ )
  3183. -        tmp1_re->ve[j] = 0.0;
  3184. -        j = limit;
  3185. -        while ( j >= 0 )
  3186. -        {
  3187. -        if ( j > 0 && T->me[j][j-1] != 0.0 )
  3188. -        {   /* 2 x 2 diagonal block */
  3189. -            /* printf("checkpoint A\n"); */
  3190. -            val1_re = tmp1_re->ve[j-1] -
  3191. -              __ip__(&(tmp1_re->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  3192. -            /* printf("checkpoint B\n"); */
  3193. -            val1_im = tmp1_im->ve[j-1] -
  3194. -              __ip__(&(tmp1_im->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  3195. -            /* printf("checkpoint C\n"); */
  3196. -            val2_re = tmp1_re->ve[j] -
  3197. -              __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  3198. -            /* printf("checkpoint D\n"); */
  3199. -            val2_im = tmp1_im->ve[j] -
  3200. -              __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  3201. -            /* printf("checkpoint E\n"); */
  3202. -            
  3203. -            t11_re = T_me[j-1][j-1] - l_re;
  3204. -            t11_im = - l_im;
  3205. -            t22_re = T_me[j][j] - l_re;
  3206. -            t22_im = - l_im;
  3207. -            t12 = T_me[j-1][j];
  3208. -            t21 = T_me[j][j-1];
  3209. -
  3210. -            scale =  fabs(T_me[j-1][j-1]) + fabs(T_me[j][j]) +
  3211. -            fabs(t12) + fabs(t21) + fabs(l_re) + fabs(l_im);
  3212. -
  3213. -            det_re = t11_re*t22_re - t11_im*t22_im - t12*t21;
  3214. -            det_im = t11_re*t22_im + t11_im*t22_re;
  3215. -            magdet = det_re*det_re+det_im*det_im;
  3216. -            if ( sqrt(magdet) < MACHEPS*scale )
  3217. -            {
  3218. -                det_re = MACHEPS*scale;
  3219. -            magdet = det_re*det_re+det_im*det_im;
  3220. -            }
  3221. -            invdet_re =   det_re/magdet;
  3222. -            invdet_im = - det_im/magdet;
  3223. -            tmp_val1_re = t22_re*val1_re-t22_im*val1_im-t12*val2_re;
  3224. -            tmp_val1_im = t22_im*val1_re+t22_re*val1_im-t12*val2_im;
  3225. -            tmp_val2_re = t11_re*val2_re-t11_im*val2_im-t21*val1_re;
  3226. -            tmp_val2_im = t11_im*val2_re+t11_re*val2_im-t21*val1_im;
  3227. -            tmp1_re->ve[j-1] = invdet_re*tmp_val1_re -
  3228. -                    invdet_im*tmp_val1_im;
  3229. -            tmp1_im->ve[j-1] = invdet_im*tmp_val1_re +
  3230. -                    invdet_re*tmp_val1_im;
  3231. -            tmp1_re->ve[j]   = invdet_re*tmp_val2_re -
  3232. -                    invdet_im*tmp_val2_im;
  3233. -            tmp1_im->ve[j]   = invdet_im*tmp_val2_re +
  3234. -                    invdet_re*tmp_val2_im;
  3235. -            j -= 2;
  3236. -            }
  3237. -            else
  3238. -        {
  3239. -            t11_re = T_me[j][j] - l_re;
  3240. -            t11_im = - l_im;
  3241. -            magdet = t11_re*t11_re + t11_im*t11_im;
  3242. -            scale = fabs(T_me[j][j]) + fabs(l_re);
  3243. -            if ( sqrt(magdet) < MACHEPS*scale )
  3244. -            {
  3245. -                t11_re = MACHEPS*scale;
  3246. -            magdet = t11_re*t11_re + t11_im*t11_im;
  3247. -            }
  3248. -            invdet_re =   t11_re/magdet;
  3249. -            invdet_im = - t11_im/magdet;
  3250. -            /* printf("checkpoint F\n"); */
  3251. -            val1_re = tmp1_re->ve[j] -
  3252. -              __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  3253. -            /* printf("checkpoint G\n"); */
  3254. -            val1_im = tmp1_im->ve[j] -
  3255. -              __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  3256. -            /* printf("checkpoint H\n"); */
  3257. -            tmp1_re->ve[j] = invdet_re*val1_re - invdet_im*val1_im;
  3258. -            tmp1_im->ve[j] = invdet_im*val1_re + invdet_re*val1_im;
  3259. -            j -= 1;
  3260. -        }
  3261. -        }
  3262. -
  3263. -        norm = v_norm_inf(tmp1_re) + v_norm_inf(tmp1_im);
  3264. -        sv_mlt(1/norm,tmp1_re,tmp1_re);
  3265. -        if ( l_im != 0.0 )
  3266. -        sv_mlt(1/norm,tmp1_im,tmp1_im);
  3267. -        mv_mlt(Q,tmp1_re,tmp2_re);
  3268. -        if ( l_im != 0.0 )
  3269. -        mv_mlt(Q,tmp1_im,tmp2_im);
  3270. -        if ( l_im != 0.0 )
  3271. -        norm = sqrt(in_prod(tmp2_re,tmp2_re)+in_prod(tmp2_im,tmp2_im));
  3272. -        else
  3273. -        norm = v_norm2(tmp2_re);
  3274. -        sv_mlt(1/norm,tmp2_re,tmp2_re);
  3275. -        if ( l_im != 0.0 )
  3276. -        sv_mlt(1/norm,tmp2_im,tmp2_im);
  3277. -
  3278. -        if ( l_im != 0.0 )
  3279. -        {
  3280. -        if ( ! X_im )
  3281. -        error(E_NULL,"schur_vecs");
  3282. -        set_col(X_re,i,tmp2_re);
  3283. -        set_col(X_im,i,tmp2_im);
  3284. -        sv_mlt(-1.0,tmp2_im,tmp2_im);
  3285. -        set_col(X_re,i+1,tmp2_re);
  3286. -        set_col(X_im,i+1,tmp2_im);
  3287. -        i += 2;
  3288. -        }
  3289. -        else
  3290. -        {
  3291. -        set_col(X_re,i,tmp2_re);
  3292. -        if ( X_im != MNULL )
  3293. -            set_col(X_im,i,tmp1_im);    /* zero vector */
  3294. -        i += 1;
  3295. -        }
  3296. -    }
  3297. -
  3298. -    return X_re;
  3299. -}
  3300. -
  3301. //GO.SYSIN DD schur.c
  3302. echo svd.c 1>&2
  3303. sed >svd.c <<'//GO.SYSIN DD svd.c' 's/^-//'
  3304. -
  3305. -/**************************************************************************
  3306. -**
  3307. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3308. -**
  3309. -**                 Meschach Library
  3310. -** 
  3311. -** This Meschach Library is provided "as is" without any express 
  3312. -** or implied warranty of any kind with respect to this software. 
  3313. -** In particular the authors shall not be liable for any direct, 
  3314. -** indirect, special, incidental or consequential damages arising 
  3315. -** in any way from use of the software.
  3316. -** 
  3317. -** Everyone is granted permission to copy, modify and redistribute this
  3318. -** Meschach Library, provided:
  3319. -**  1.  All copies contain this copyright notice.
  3320. -**  2.  All modified copies shall carry a notice stating who
  3321. -**      made the last modification and the date of such modification.
  3322. -**  3.  No charge is made for this software or works derived from it.  
  3323. -**      This clause shall not be construed as constraining other software
  3324. -**      distributed on the same medium as this software, nor is a
  3325. -**      distribution fee considered a charge.
  3326. -**
  3327. -***************************************************************************/
  3328. -
  3329. -
  3330. -/*
  3331. -    File containing routines for computing the SVD of matrices
  3332. -*/
  3333. -
  3334. -#include    <stdio.h>
  3335. -#include    <math.h>
  3336. -#include    "matrix.h"
  3337. -#include        "matrix2.h"
  3338. -
  3339. -
  3340. -static char rcsid[] = "$Id: svd.c,v 1.6 1994/01/13 05:44:16 des Exp $";
  3341. -
  3342. -
  3343. -
  3344. -#define    sgn(x)    ((x) >= 0 ? 1 : -1)
  3345. -#define    MAX_STACK    100
  3346. -
  3347. -/* fixsvd -- fix minor details about SVD
  3348. -    -- make singular values non-negative
  3349. -    -- sort singular values in decreasing order
  3350. -    -- variables as for bisvd()
  3351. -    -- no argument checking */
  3352. -static void    fixsvd(d,U,V)
  3353. -VEC    *d;
  3354. -MAT    *U, *V;
  3355. -{
  3356. -    int        i, j, k, l, r, stack[MAX_STACK], sp;
  3357. -    Real    tmp, v;
  3358. -
  3359. -    /* make singular values non-negative */
  3360. -    for ( i = 0; i < d->dim; i++ )
  3361. -    if ( d->ve[i] < 0.0 )
  3362. -    {
  3363. -        d->ve[i] = - d->ve[i];
  3364. -        if ( U != MNULL )
  3365. -        for ( j = 0; j < U->m; j++ )
  3366. -            U->me[i][j] = - U->me[i][j];
  3367. -    }
  3368. -
  3369. -    /* sort singular values */
  3370. -    /* nonrecursive implementation of quicksort due to R.Sedgewick,
  3371. -       "Algorithms in C", p. 122 (1990) */
  3372. -    sp = -1;
  3373. -    l = 0;    r = d->dim - 1;
  3374. -    for ( ; ; )
  3375. -    {
  3376. -    while ( r > l )
  3377. -    {
  3378. -        /* i = partition(d->ve,l,r) */
  3379. -        v = d->ve[r];
  3380. -
  3381. -        i = l - 1;        j = r;
  3382. -        for ( ; ; )
  3383. -        {    /* inequalities are "backwards" for **decreasing** order */
  3384. -        while ( d->ve[++i] > v )
  3385. -            ;
  3386. -        while ( d->ve[--j] < v )
  3387. -            ;
  3388. -        if ( i >= j )
  3389. -            break;
  3390. -        /* swap entries in d->ve */
  3391. -        tmp = d->ve[i];      d->ve[i] = d->ve[j];    d->ve[j] = tmp;
  3392. -        /* swap rows of U & V as well */
  3393. -        if ( U != MNULL )
  3394. -            for ( k = 0; k < U->n; k++ )
  3395. -            {
  3396. -            tmp = U->me[i][k];
  3397. -            U->me[i][k] = U->me[j][k];
  3398. -            U->me[j][k] = tmp;
  3399. -            }
  3400. -        if ( V != MNULL )
  3401. -            for ( k = 0; k < V->n; k++ )
  3402. -            {
  3403. -            tmp = V->me[i][k];
  3404. -            V->me[i][k] = V->me[j][k];
  3405. -            V->me[j][k] = tmp;
  3406. -            }
  3407. -        }
  3408. -        tmp = d->ve[i];    d->ve[i] = d->ve[r];    d->ve[r] = tmp;
  3409. -        if ( U != MNULL )
  3410. -        for ( k = 0; k < U->n; k++ )
  3411. -        {
  3412. -            tmp = U->me[i][k];
  3413. -            U->me[i][k] = U->me[r][k];
  3414. -            U->me[r][k] = tmp;
  3415. -        }
  3416. -        if ( V != MNULL )
  3417. -        for ( k = 0; k < V->n; k++ )
  3418. -        {
  3419. -            tmp = V->me[i][k];
  3420. -            V->me[i][k] = V->me[r][k];
  3421. -            V->me[r][k] = tmp;
  3422. -        }
  3423. -        /* end i = partition(...) */
  3424. -        if ( i - l > r - i )
  3425. -        {    stack[++sp] = l;    stack[++sp] = i-1;    l = i+1;    }
  3426. -        else
  3427. -        {    stack[++sp] = i+1;  stack[++sp] = r;    r = i-1;    }
  3428. -    }
  3429. -    if ( sp < 0 )
  3430. -        break;
  3431. -    r = stack[sp--];    l = stack[sp--];
  3432. -    }
  3433. -}
  3434. -
  3435. -
  3436. -/* bisvd -- svd of a bidiagonal m x n matrix represented by d (diagonal) and
  3437. -            f (super-diagonals)
  3438. -    -- returns with d set to the singular values, f zeroed
  3439. -    -- if U, V non-NULL, the orthogonal operations are accumulated
  3440. -        in U, V; if U, V == I on entry, then SVD == U^T.A.V
  3441. -        where A is initial matrix
  3442. -    -- returns d on exit */
  3443. -VEC    *bisvd(d,f,U,V)
  3444. -VEC    *d, *f;
  3445. -MAT    *U, *V;
  3446. -{
  3447. -    int    i, j, n;
  3448. -    int    i_min, i_max, split;
  3449. -    Real    c, s, shift, size, z;
  3450. -    Real    d_tmp, diff, t11, t12, t22, *d_ve, *f_ve;
  3451. -
  3452. -    if ( ! d || ! f )
  3453. -        error(E_NULL,"bisvd");
  3454. -    if ( d->dim != f->dim + 1 )
  3455. -        error(E_SIZES,"bisvd");
  3456. -    n = d->dim;
  3457. -    if ( ( U && U->n < n ) || ( V && V->m < n ) )
  3458. -        error(E_SIZES,"bisvd");
  3459. -    if ( ( U && U->m != U->n ) || ( V && V->m != V->n ) )
  3460. -        error(E_SQUARE,"bisvd");
  3461. -
  3462. -
  3463. -    if ( n == 1 )
  3464. -        return d;
  3465. -    d_ve = d->ve;    f_ve = f->ve;
  3466. -
  3467. -    size = v_norm_inf(d) + v_norm_inf(f);
  3468. -
  3469. -    i_min = 0;
  3470. -    while ( i_min < n )    /* outer while loop */
  3471. -    {
  3472. -        /* find i_max to suit;
  3473. -        submatrix i_min..i_max should be irreducible */
  3474. -        i_max = n - 1;
  3475. -        for ( i = i_min; i < n - 1; i++ )
  3476. -        if ( d_ve[i] == 0.0 || f_ve[i] == 0.0 )
  3477. -        {   i_max = i;
  3478. -            if ( f_ve[i] != 0.0 )
  3479. -            {
  3480. -            /* have to ``chase'' f[i] element out of matrix */
  3481. -            z = f_ve[i];    f_ve[i] = 0.0;
  3482. -            for ( j = i; j < n-1 && z != 0.0; j++ )
  3483. -            {
  3484. -                givens(d_ve[j+1],z, &c, &s);
  3485. -                s = -s;
  3486. -                d_ve[j+1] =  c*d_ve[j+1] - s*z;
  3487. -                if ( j+1 < n-1 )
  3488. -                {
  3489. -                z         = s*f_ve[j+1];
  3490. -                f_ve[j+1] = c*f_ve[j+1];
  3491. -                }
  3492. -                if ( U )
  3493. -                rot_rows(U,i,j+1,c,s,U);
  3494. -            }
  3495. -            }
  3496. -            break;
  3497. -        }
  3498. -        if ( i_max <= i_min )
  3499. -        {
  3500. -        i_min = i_max + 1;
  3501. -        continue;
  3502. -        }
  3503. -        /* printf("bisvd: i_min = %d, i_max = %d\n",i_min,i_max); */
  3504. -
  3505. -        split = FALSE;
  3506. -        while ( ! split )
  3507. -        {
  3508. -        /* compute shift */
  3509. -        t11 = d_ve[i_max-1]*d_ve[i_max-1] +
  3510. -            (i_max > i_min+1 ? f_ve[i_max-2]*f_ve[i_max-2] : 0.0);
  3511. -        t12 = d_ve[i_max-1]*f_ve[i_max-1];
  3512. -        t22 = d_ve[i_max]*d_ve[i_max] + f_ve[i_max-1]*f_ve[i_max-1];
  3513. -        /* use e-val of [[t11,t12],[t12,t22]] matrix
  3514. -                closest to t22 */
  3515. -        diff = (t11-t22)/2;
  3516. -        shift = t22 - t12*t12/(diff +
  3517. -            sgn(diff)*sqrt(diff*diff+t12*t12));
  3518. -
  3519. -        /* initial Givens' rotation */
  3520. -        givens(d_ve[i_min]*d_ve[i_min]-shift,
  3521. -            d_ve[i_min]*f_ve[i_min], &c, &s);
  3522. -
  3523. -        /* do initial Givens' rotations */
  3524. -        d_tmp         = c*d_ve[i_min] + s*f_ve[i_min];
  3525. -        f_ve[i_min]   = c*f_ve[i_min] - s*d_ve[i_min];
  3526. -        d_ve[i_min]   = d_tmp;
  3527. -        z             = s*d_ve[i_min+1];
  3528. -        d_ve[i_min+1] = c*d_ve[i_min+1];
  3529. -        if ( V )
  3530. -            rot_rows(V,i_min,i_min+1,c,s,V);
  3531. -        /* 2nd Givens' rotation */
  3532. -        givens(d_ve[i_min],z, &c, &s);
  3533. -        d_ve[i_min]   = c*d_ve[i_min] + s*z;
  3534. -        d_tmp         = c*d_ve[i_min+1] - s*f_ve[i_min];
  3535. -        f_ve[i_min]   = s*d_ve[i_min+1] + c*f_ve[i_min];
  3536. -        d_ve[i_min+1] = d_tmp;
  3537. -        if ( i_min+1 < i_max )
  3538. -        {
  3539. -            z             = s*f_ve[i_min+1];
  3540. -            f_ve[i_min+1] = c*f_ve[i_min+1];
  3541. -        }
  3542. -        if ( U )
  3543. -            rot_rows(U,i_min,i_min+1,c,s,U);
  3544. -
  3545. -        for ( i = i_min+1; i < i_max; i++ )
  3546. -        {
  3547. -            /* get Givens' rotation for zeroing z */
  3548. -            givens(f_ve[i-1],z, &c, &s);
  3549. -            f_ve[i-1] = c*f_ve[i-1] + s*z;
  3550. -            d_tmp     = c*d_ve[i] + s*f_ve[i];
  3551. -            f_ve[i]   = c*f_ve[i] - s*d_ve[i];
  3552. -            d_ve[i]   = d_tmp;
  3553. -            z         = s*d_ve[i+1];
  3554. -            d_ve[i+1] = c*d_ve[i+1];
  3555. -            if ( V )
  3556. -            rot_rows(V,i,i+1,c,s,V);
  3557. -            /* get 2nd Givens' rotation */
  3558. -            givens(d_ve[i],z, &c, &s);
  3559. -            d_ve[i]   = c*d_ve[i] + s*z;
  3560. -            d_tmp     = c*d_ve[i+1] - s*f_ve[i];
  3561. -            f_ve[i]   = c*f_ve[i] + s*d_ve[i+1];
  3562. -            d_ve[i+1] = d_tmp;
  3563. -            if ( i+1 < i_max )
  3564. -            {
  3565. -            z         = s*f_ve[i+1];
  3566. -            f_ve[i+1] = c*f_ve[i+1];
  3567. -            }
  3568. -            if ( U )
  3569. -            rot_rows(U,i,i+1,c,s,U);
  3570. -        }
  3571. -        /* should matrix be split? */
  3572. -        for ( i = i_min; i < i_max; i++ )
  3573. -            if ( fabs(f_ve[i]) <
  3574. -                MACHEPS*(fabs(d_ve[i])+fabs(d_ve[i+1])) )
  3575. -            {
  3576. -            split = TRUE;
  3577. -            f_ve[i] = 0.0;
  3578. -            }
  3579. -            else if ( fabs(d_ve[i]) < MACHEPS*size )
  3580. -            {
  3581. -            split = TRUE;
  3582. -            d_ve[i] = 0.0;
  3583. -            }
  3584. -            /* printf("bisvd: d =\n");    v_output(d); */
  3585. -            /* printf("bisvd: f = \n");    v_output(f); */
  3586. -        }
  3587. -    }
  3588. -    fixsvd(d,U,V);
  3589. -
  3590. -    return d;
  3591. -}
  3592. -
  3593. -/* bifactor -- perform preliminary factorisation for bisvd
  3594. -    -- updates U and/or V, which ever is not NULL */
  3595. -MAT    *bifactor(A,U,V)
  3596. -MAT    *A, *U, *V;
  3597. -{
  3598. -    int    k;
  3599. -    static VEC    *tmp1=VNULL, *tmp2=VNULL;
  3600. -    Real    beta;
  3601. -
  3602. -    if ( ! A )
  3603. -        error(E_NULL,"bifactor");
  3604. -    if ( ( U && ( U->m != U->n ) ) || ( V && ( V->m != V->n ) ) )
  3605. -        error(E_SQUARE,"bifactor");
  3606. -    if ( ( U && U->m != A->m ) || ( V && V->m != A->n ) )
  3607. -        error(E_SIZES,"bifactor");
  3608. -    tmp1 = v_resize(tmp1,A->m);
  3609. -    tmp2 = v_resize(tmp2,A->n);
  3610. -    MEM_STAT_REG(tmp1,TYPE_VEC);
  3611. -    MEM_STAT_REG(tmp2,TYPE_VEC);
  3612. -
  3613. -    if ( A->m >= A->n )
  3614. -        for ( k = 0; k < A->n; k++ )
  3615. -        {
  3616. -        get_col(A,k,tmp1);
  3617. -        hhvec(tmp1,k,&beta,tmp1,&(A->me[k][k]));
  3618. -        hhtrcols(A,k,k+1,tmp1,beta);
  3619. -        if ( U )
  3620. -            hhtrcols(U,k,0,tmp1,beta);
  3621. -        if ( k+1 >= A->n )
  3622. -            continue;
  3623. -        get_row(A,k,tmp2);
  3624. -        hhvec(tmp2,k+1,&beta,tmp2,&(A->me[k][k+1]));
  3625. -        hhtrrows(A,k+1,k+1,tmp2,beta);
  3626. -        if ( V )
  3627. -            hhtrcols(V,k+1,0,tmp2,beta);
  3628. -        }
  3629. -    else
  3630. -        for ( k = 0; k < A->m; k++ )
  3631. -        {
  3632. -        get_row(A,k,tmp2);
  3633. -        hhvec(tmp2,k,&beta,tmp2,&(A->me[k][k]));
  3634. -        hhtrrows(A,k+1,k,tmp2,beta);
  3635. -        if ( V )
  3636. -            hhtrcols(V,k,0,tmp2,beta);
  3637. -        if ( k+1 >= A->m )
  3638. -            continue;
  3639. -        get_col(A,k,tmp1);
  3640. -        hhvec(tmp1,k+1,&beta,tmp1,&(A->me[k+1][k]));
  3641. -        hhtrcols(A,k+1,k+1,tmp1,beta);
  3642. -        if ( U )
  3643. -            hhtrcols(U,k+1,0,tmp1,beta);
  3644. -        }
  3645. -
  3646. -    return A;
  3647. -}
  3648. -
  3649. -/* svd -- returns vector of singular values in d
  3650. -    -- also updates U and/or V, if one or the other is non-NULL
  3651. -    -- destroys A */
  3652. -VEC    *svd(A,U,V,d)
  3653. -MAT    *A, *U, *V;
  3654. -VEC    *d;
  3655. -{
  3656. -    static VEC    *f=VNULL;
  3657. -    int    i, limit;
  3658. -    MAT    *A_tmp;
  3659. -
  3660. -    if ( ! A )
  3661. -        error(E_NULL,"svd");
  3662. -    if ( ( U && ( U->m != U->n ) ) || ( V && ( V->m != V->n ) ) )
  3663. -        error(E_SQUARE,"svd");
  3664. -    if ( ( U && U->m != A->m ) || ( V && V->m != A->n ) )
  3665. -        error(E_SIZES,"svd");
  3666. -
  3667. -    A_tmp = m_copy(A,MNULL);
  3668. -    if ( U != MNULL )
  3669. -        m_ident(U);
  3670. -    if ( V != MNULL )
  3671. -        m_ident(V);
  3672. -    limit = min(A_tmp->m,A_tmp->n);
  3673. -    d = v_resize(d,limit);
  3674. -    f = v_resize(f,limit-1);
  3675. -    MEM_STAT_REG(f,TYPE_VEC);
  3676. -
  3677. -    bifactor(A_tmp,U,V);
  3678. -    if ( A_tmp->m >= A_tmp->n )
  3679. -        for ( i = 0; i < limit; i++ )
  3680. -        {
  3681. -        d->ve[i] = A_tmp->me[i][i];
  3682. -        if ( i+1 < limit )
  3683. -            f->ve[i] = A_tmp->me[i][i+1];
  3684. -        }
  3685. -    else
  3686. -        for ( i = 0; i < limit; i++ )
  3687. -        {
  3688. -        d->ve[i] = A_tmp->me[i][i];
  3689. -        if ( i+1 < limit )
  3690. -            f->ve[i] = A_tmp->me[i+1][i];
  3691. -        }
  3692. -
  3693. -
  3694. -    if ( A_tmp->m >= A_tmp->n )
  3695. -        bisvd(d,f,U,V);
  3696. -    else
  3697. -        bisvd(d,f,V,U);
  3698. -
  3699. -    M_FREE(A_tmp);
  3700. -
  3701. -    return d;
  3702. -}
  3703. -
  3704. //GO.SYSIN DD svd.c
  3705. echo fft.c 1>&2
  3706. sed >fft.c <<'//GO.SYSIN DD fft.c' 's/^-//'
  3707. -
  3708. -/**************************************************************************
  3709. -**
  3710. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3711. -**
  3712. -**                 Meschach Library
  3713. -** 
  3714. -** This Meschach Library is provided "as is" without any express 
  3715. -** or implied warranty of any kind with respect to this software. 
  3716. -** In particular the authors shall not be liable for any direct, 
  3717. -** indirect, special, incidental or consequential damages arising 
  3718. -** in any way from use of the software.
  3719. -** 
  3720. -** Everyone is granted permission to copy, modify and redistribute this
  3721. -** Meschach Library, provided:
  3722. -**  1.  All copies contain this copyright notice.
  3723. -**  2.  All modified copies shall carry a notice stating who
  3724. -**      made the last modification and the date of such modification.
  3725. -**  3.  No charge is made for this software or works derived from it.  
  3726. -**      This clause shall not be construed as constraining other software
  3727. -**      distributed on the same medium as this software, nor is a
  3728. -**      distribution fee considered a charge.
  3729. -**
  3730. -***************************************************************************/
  3731. -
  3732. -
  3733. -/*
  3734. -    Fast Fourier Transform routine
  3735. -    Loosely based on the Fortran routine in Rabiner & Gold's
  3736. -    "Digital Signal Processing"
  3737. -*/
  3738. -
  3739. -static char rcsid[] = "$Id: fft.c,v 1.3 1994/01/13 05:45:33 des Exp $";
  3740. -
  3741. -#include        <stdio.h>
  3742. -#include        <math.h>
  3743. -#include        "matrix.h"
  3744. -#include        "matrix2.h"
  3745. -
  3746. -
  3747. -/* fft -- d.i.t. fast Fourier transform 
  3748. -        -- radix-2 FFT only
  3749. -        -- vector extended to a power of 2 */
  3750. -void    fft(x_re,x_im)
  3751. -VEC     *x_re, *x_im;
  3752. -{
  3753. -    int         i, ip, j, k, li, n, length;
  3754. -    Real      *xr, *xi;
  3755. -    Real    theta, pi = 3.1415926535897932384;
  3756. -    Real      w_re, w_im, u_re, u_im, t_re, t_im;
  3757. -    Real      tmp, tmpr, tmpi;
  3758. -
  3759. -    if ( ! x_re || ! x_im )
  3760. -        error(E_NULL,"fft");
  3761. -    if ( x_re->dim != x_im->dim )
  3762. -        error(E_SIZES,"fft");
  3763. -
  3764. -    n = 1;
  3765. -    while ( x_re->dim > n )
  3766. -        n *= 2;
  3767. -    x_re = v_resize(x_re,n);
  3768. -    x_im = v_resize(x_im,n);
  3769. -    printf("# fft: x_re =\n");  v_output(x_re);
  3770. -    printf("# fft: x_im =\n");  v_output(x_im);
  3771. -    xr   = x_re->ve;
  3772. -    xi   = x_im->ve;
  3773. -
  3774. -    /* Decimation in time (DIT) algorithm */
  3775. -    j = 0;
  3776. -    for ( i = 0; i < n-1; i++ )
  3777. -    {
  3778. -        if ( i < j )
  3779. -        {
  3780. -            tmp = xr[i];
  3781. -            xr[i] = xr[j];
  3782. -            xr[j] = tmp;
  3783. -            tmp = xi[i];
  3784. -            xi[i] = xi[j];
  3785. -            xi[j] = tmp;
  3786. -        }
  3787. -        k = n / 2;
  3788. -        while ( k <= j )
  3789. -        {
  3790. -            j -= k;
  3791. -            k /= 2;
  3792. -        }
  3793. -        j += k;
  3794. -    }
  3795. -
  3796. -    /* Actual FFT */
  3797. -    for ( li = 1; li < n; li *= 2 )
  3798. -    {
  3799. -        length = 2*li;
  3800. -        theta  = pi/li;
  3801. -        u_re   = 1.0;
  3802. -        u_im   = 0.0;
  3803. -        if ( li == 1 )
  3804. -        {
  3805. -            w_re = -1.0;
  3806. -            w_im =  0.0;
  3807. -        }
  3808. -        else if ( li == 2 )
  3809. -        {
  3810. -            w_re =  0.0;
  3811. -            w_im =  1.0;
  3812. -        }
  3813. -        else
  3814. -        {
  3815. -            w_re = cos(theta);
  3816. -            w_im = sin(theta);
  3817. -        }
  3818. -        for ( j = 0; j < li; j++ )
  3819. -        {
  3820. -            for ( i =  j; i < n; i += length )
  3821. -            {
  3822. -                ip = i + li;
  3823. -                /* step 1 */
  3824. -                t_re = xr[ip]*u_re - xi[ip]*u_im;
  3825. -                t_im = xr[ip]*u_im + xi[ip]*u_re;
  3826. -                /* step 2 */
  3827. -                xr[ip] = xr[i]  - t_re;
  3828. -                xi[ip] = xi[i]  - t_im;
  3829. -                /* step 3 */
  3830. -                xr[i] += t_re;
  3831. -                xi[i] += t_im;
  3832. -            }
  3833. -            tmpr = u_re*w_re - u_im*w_im;
  3834. -            tmpi = u_im*w_re + u_re*w_im;
  3835. -            u_re = tmpr;
  3836. -            u_im = tmpi;
  3837. -        }
  3838. -    }
  3839. -}
  3840. -
  3841. -/* ifft -- inverse FFT using the same interface as fft() */
  3842. -void    ifft(x_re,x_im)
  3843. -VEC    *x_re, *x_im;
  3844. -{
  3845. -    /* we just use complex conjugates */
  3846. -
  3847. -    sv_mlt(-1.0,x_im,x_im);
  3848. -    fft(x_re,x_im);
  3849. -    sv_mlt(-1.0/((double)(x_re->dim)),x_im,x_im);
  3850. -}
  3851. //GO.SYSIN DD fft.c
  3852. echo mfunc.c 1>&2
  3853. sed >mfunc.c <<'//GO.SYSIN DD mfunc.c' 's/^-//'
  3854. -
  3855. -/**************************************************************************
  3856. -**
  3857. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3858. -**
  3859. -**                 Meschach Library
  3860. -** 
  3861. -** This Meschach Library is provided "as is" without any express 
  3862. -** or implied warranty of any kind with respect to this software. 
  3863. -** In particular the authors shall not be liable for any direct, 
  3864. -** indirect, special, incidental or consequential damages arising 
  3865. -** in any way from use of the software.
  3866. -** 
  3867. -** Everyone is granted permission to copy, modify and redistribute this
  3868. -** Meschach Library, provided:
  3869. -**  1.  All copies contain this copyright notice.
  3870. -**  2.  All modified copies shall carry a notice stating who
  3871. -**      made the last modification and the date of such modification.
  3872. -**  3.  No charge is made for this software or works derived from it.  
  3873. -**      This clause shall not be construed as constraining other software
  3874. -**      distributed on the same medium as this software, nor is a
  3875. -**      distribution fee considered a charge.
  3876. -**
  3877. -***************************************************************************/
  3878. -
  3879. -
  3880. -/*
  3881. -  This file contains routines for computing functions of matrices
  3882. -  especially polynomials and exponential functions
  3883. -  Copyright (C) Teresa Leyk and David Stewart, 1993
  3884. -  */
  3885. -
  3886. -#include <stdio.h>
  3887. -#include <math.h>
  3888. -#include "matrix.h"
  3889. -#include "matrix2.h"
  3890. -
  3891. -static char    rcsid[] = "$Id: mfunc.c,v 1.1 1994/01/13 05:33:58 des Exp $";
  3892. -
  3893. -
  3894. -
  3895. -/* _m_pow -- computes integer powers of a square matrix A, A^p
  3896. -   -- uses tmp as temporary workspace */
  3897. -MAT    *_m_pow(A, p, tmp, out)
  3898. -MAT    *A, *tmp, *out;
  3899. -int    p;
  3900. -{
  3901. -   int        it_cnt, k, max_bit;
  3902. -   
  3903. -   /*
  3904. -     File containing routines for evaluating matrix functions
  3905. -     esp. the exponential function
  3906. -     */
  3907. -
  3908. -#define    Z(k)    (((k) & 1) ? tmp : out)
  3909. -   
  3910. -   if ( ! A )
  3911. -     error(E_NULL,"_m_pow");
  3912. -   if ( A->m != A->n )
  3913. -     error(E_SQUARE,"_m_pow");
  3914. -   if ( p < 0 )
  3915. -     error(E_NEG,"_m_pow");
  3916. -   out = m_resize(out,A->m,A->n);
  3917. -   tmp = m_resize(tmp,A->m,A->n);
  3918. -   
  3919. -   if ( p == 0 )
  3920. -     m_ident(out);
  3921. -   else if ( p > 0 )
  3922. -   {
  3923. -      it_cnt = 1;
  3924. -      for ( max_bit = 0; ; max_bit++ )
  3925. -    if ( (p >> (max_bit+1)) == 0 )
  3926. -      break;
  3927. -      tmp = m_copy(A,tmp);
  3928. -      
  3929. -      for ( k = 0; k < max_bit; k++ )
  3930. -      {
  3931. -     m_mlt(Z(it_cnt),Z(it_cnt),Z(it_cnt+1));
  3932. -     it_cnt++;
  3933. -     if ( p & (1 << (max_bit-1)) )
  3934. -     {
  3935. -        m_mlt(A,Z(it_cnt),Z(it_cnt+1));
  3936. -        /* m_copy(Z(it_cnt),out); */
  3937. -        it_cnt++;
  3938. -     }
  3939. -     p <<= 1;
  3940. -      }
  3941. -      if (it_cnt & 1)
  3942. -    out = m_copy(Z(it_cnt),out);
  3943. -   }
  3944. -
  3945. -   return out;
  3946. -
  3947. -#undef Z   
  3948. -}
  3949. -
  3950. -/* m_pow -- computes integer powers of a square matrix A, A^p */
  3951. -MAT    *m_pow(A, p, out)
  3952. -MAT    *A, *out;
  3953. -int    p;
  3954. -{
  3955. -   static MAT    *wkspace = MNULL;
  3956. -   
  3957. -   if ( ! A )
  3958. -     error(E_NULL,"m_pow");
  3959. -   if ( A->m != A->n )
  3960. -     error(E_SQUARE,"m_pow");
  3961. -   
  3962. -   wkspace = m_resize(wkspace,A->m,A->n);
  3963. -   MEM_STAT_REG(wkspace,TYPE_MAT);
  3964. -   
  3965. -   return _m_pow(A, p, wkspace, out);
  3966. -   
  3967. -}
  3968. -
  3969. -/**************************************************/
  3970. -
  3971. -/* _m_exp -- compute matrix exponential of A and save it in out
  3972. -   -- uses Pade approximation followed by repeated squaring
  3973. -   -- eps is the tolerance used for the Pade approximation 
  3974. -   -- A is not changed
  3975. -   -- q_out - degree of the Pade approximation (q_out,q_out)
  3976. -   -- j_out - the power of 2 for scaling the matrix A
  3977. -              such that ||A/2^j_out|| <= 0.5
  3978. -*/
  3979. -MAT *_m_exp(A,eps,out,q_out,j_out)
  3980. -MAT *A,*out;
  3981. -double eps;
  3982. -int *q_out, *j_out;
  3983. -{
  3984. -   static MAT *D = MNULL, *Apow = MNULL, *N = MNULL, *Y = MNULL;
  3985. -   static VEC *c1 = VNULL, *tmp = VNULL;
  3986. -   VEC y0, y1;  /* additional structures */
  3987. -   static PERM *pivot = PNULL;
  3988. -   int j, k, l, q, r, s, j2max, t;
  3989. -   double inf_norm, eqq, power2, c, sign;
  3990. -   
  3991. -   if ( ! A )
  3992. -     error(E_SIZES,"_m_exp");
  3993. -   if ( A->m != A->n )
  3994. -     error(E_SIZES,"_m_exp");
  3995. -   if ( A == out )
  3996. -     error(E_INSITU,"_m_exp");
  3997. -   if ( eps < 0.0 )
  3998. -     error(E_RANGE,"_m_exp");
  3999. -   else if (eps == 0.0)
  4000. -     eps = MACHEPS;
  4001. -      
  4002. -   N = m_resize(N,A->m,A->n);
  4003. -   D = m_resize(D,A->m,A->n);
  4004. -   Apow = m_resize(Apow,A->m,A->n);
  4005. -   out = m_resize(out,A->m,A->n);
  4006. -
  4007. -   MEM_STAT_REG(N,TYPE_MAT);
  4008. -   MEM_STAT_REG(D,TYPE_MAT);
  4009. -   MEM_STAT_REG(Apow,TYPE_MAT);
  4010. -   
  4011. -   /* normalise A to have ||A||_inf <= 1 */
  4012. -   inf_norm = m_norm_inf(A);
  4013. -   if (inf_norm <= 0.0) {
  4014. -      m_ident(out);
  4015. -      *q_out = -1;
  4016. -      *j_out = 0;
  4017. -      return out;
  4018. -   }
  4019. -   else {
  4020. -      j2max = floor(1+log(inf_norm)/log(2.0));
  4021. -      j2max = max(0, j2max);
  4022. -   }
  4023. -   
  4024. -   power2 = 1.0;
  4025. -   for ( k = 1; k <= j2max; k++ )
  4026. -     power2 *= 2;
  4027. -   power2 = 1.0/power2;
  4028. -   if ( j2max > 0 )
  4029. -     sm_mlt(power2,A,A);
  4030. -   
  4031. -   /* compute order for polynomial approximation */
  4032. -   eqq = 1.0/6.0;
  4033. -   for ( q = 1; eqq > eps; q++ )
  4034. -     eqq /= 16.0*(2.0*q+1.0)*(2.0*q+3.0);
  4035. -   
  4036. -   /* construct vector of coefficients */
  4037. -   c1 = v_resize(c1,q+1);
  4038. -   MEM_STAT_REG(c1,TYPE_VEC);
  4039. -   c1->ve[0] = 1.0;
  4040. -   for ( k = 1; k <= q; k++ ) 
  4041. -     c1->ve[k] = c1->ve[k-1]*(q-k+1)/((2*q-k+1)*(double)k);
  4042. -   
  4043. -   tmp = v_resize(tmp,A->n);
  4044. -   MEM_STAT_REG(tmp,TYPE_VEC);
  4045. -   
  4046. -   s = (int)floor(sqrt((double)q/2.0));
  4047. -   if ( s <= 0 )  s = 1;
  4048. -   _m_pow(A,s,out,Apow);
  4049. -   r = q/s;
  4050. -   
  4051. -   Y = m_resize(Y,s,A->n);
  4052. -   MEM_STAT_REG(Y,TYPE_MAT);
  4053. -   /* y0 and y1 are pointers to rows of Y, N and D */
  4054. -   y0.dim = y0.max_dim = A->n;   
  4055. -   y1.dim = y1.max_dim = A->n;
  4056. -   
  4057. -   m_zero(Y);
  4058. -   m_zero(N);
  4059. -   m_zero(D);
  4060. -   
  4061. -   for( j = 0; j < A->n; j++ )
  4062. -   {
  4063. -      if (j > 0)
  4064. -    Y->me[0][j-1] = 0.0;
  4065. -      y0.ve = Y->me[0];
  4066. -      y0.ve[j] = 1.0;
  4067. -      for ( k = 0; k < s-1; k++ )
  4068. -      {
  4069. -     y1.ve = Y->me[k+1];
  4070. -     mv_mlt(A,&y0,&y1);
  4071. -     y0.ve = y1.ve;
  4072. -      }
  4073. -
  4074. -      y0.ve = N->me[j];
  4075. -      y1.ve = D->me[j];
  4076. -      t = s*r;
  4077. -      for ( l = 0; l <= q-t; l++ )
  4078. -      {
  4079. -     c = c1->ve[t+l];
  4080. -     sign = ((t+l) & 1) ? -1.0 : 1.0;
  4081. -     __mltadd__(y0.ve,Y->me[l],c,     Y->n);
  4082. -     __mltadd__(y1.ve,Y->me[l],c*sign,Y->n);
  4083. -      }
  4084. -      
  4085. -      for (k=1; k <= r; k++)
  4086. -      {
  4087. -     v_copy(mv_mlt(Apow,&y0,tmp),&y0);
  4088. -     v_copy(mv_mlt(Apow,&y1,tmp),&y1);
  4089. -     t = s*(r-k);
  4090. -     for (l=0; l < s; l++)
  4091. -     {
  4092. -        c = c1->ve[t+l];
  4093. -        sign = ((t+l) & 1) ? -1.0 : 1.0;
  4094. -        __mltadd__(y0.ve,Y->me[l],c,     Y->n);
  4095. -        __mltadd__(y1.ve,Y->me[l],c*sign,Y->n);
  4096. -     }
  4097. -      }
  4098. -   }
  4099. -
  4100. -   pivot = px_resize(pivot,A->m);
  4101. -   MEM_STAT_REG(pivot,TYPE_PERM);
  4102. -   
  4103. -   /* note that N and D are transposed,
  4104. -      therefore we use LUTsolve;
  4105. -      out is saved row-wise, and must be transposed 
  4106. -      after this */
  4107. -
  4108. -   LUfactor(D,pivot);
  4109. -   for (k=0; k < A->n; k++)
  4110. -   {
  4111. -      y0.ve = N->me[k];
  4112. -      y1.ve = out->me[k];
  4113. -      LUTsolve(D,pivot,&y0,&y1);
  4114. -   }
  4115. -   m_transp(out,out); 
  4116. -
  4117. -
  4118. -   /* Use recursive squaring to turn the normalised exponential to the
  4119. -      true exponential */
  4120. -
  4121. -#define Z(k)    ((k) & 1 ? Apow : out)
  4122. -
  4123. -   for( k = 1; k <= j2max; k++)
  4124. -      m_mlt(Z(k-1),Z(k-1),Z(k));
  4125. -
  4126. -   if (Z(k) == out)
  4127. -     m_copy(Apow,out);
  4128. -   
  4129. -   /* output parameters */
  4130. -   *j_out = j2max;
  4131. -   *q_out = q;
  4132. -
  4133. -   /* restore the matrix A */
  4134. -   sm_mlt(1.0/power2,A,A);
  4135. -   return out;
  4136. -
  4137. -#undef Z
  4138. -}
  4139. -
  4140. -
  4141. -/* simple interface for _m_exp */
  4142. -MAT *m_exp(A,eps,out)
  4143. -MAT *A,*out;
  4144. -double eps;
  4145. -{
  4146. -   int q_out, j_out;
  4147. -
  4148. -   return _m_exp(A,eps,out,&q_out,&j_out);
  4149. -}
  4150. -
  4151. -
  4152. -/*--------------------------------*/
  4153. -
  4154. -/* m_poly -- computes sum_i a[i].A^i, where i=0,1,...dim(a);
  4155. -   -- uses C. Van Loan's fast and memory efficient method  */
  4156. -MAT *m_poly(A,a,out)
  4157. -MAT *A,*out;
  4158. -VEC *a;
  4159. -{
  4160. -   static MAT    *Apow = MNULL, *Y = MNULL;
  4161. -   static VEC   *tmp;
  4162. -   VEC y0, y1;  /* additional vectors */
  4163. -   int j, k, l, q, r, s, t;
  4164. -   
  4165. -   if ( ! A || ! a )
  4166. -     error(E_NULL,"m_poly");
  4167. -   if ( A->m != A->n )
  4168. -     error(E_SIZES,"m_poly");
  4169. -   if ( A == out )
  4170. -     error(E_INSITU,"m_poly");
  4171. -   
  4172. -   out = m_resize(out,A->m,A->n);
  4173. -   Apow = m_resize(Apow,A->m,A->n);
  4174. -   MEM_STAT_REG(Apow,TYPE_MAT);
  4175. -   tmp = v_resize(tmp,A->n);
  4176. -   MEM_STAT_REG(tmp,TYPE_VEC);
  4177. -
  4178. -   q = a->dim - 1;
  4179. -   if ( q == 0 ) {
  4180. -      m_zero(out);
  4181. -      for (j=0; j < out->n; j++)
  4182. -    out->me[j][j] = a->ve[0];
  4183. -      return out;
  4184. -   }
  4185. -   else if ( q == 1) {
  4186. -      sm_mlt(a->ve[1],A,out);
  4187. -      for (j=0; j < out->n; j++)
  4188. -    out->me[j][j] += a->ve[0];
  4189. -      return out;
  4190. -   }
  4191. -   
  4192. -   s = (int)floor(sqrt((double)q/2.0));
  4193. -   if ( s <= 0 ) s = 1;
  4194. -   _m_pow(A,s,out,Apow);
  4195. -   r = q/s;
  4196. -   
  4197. -   Y = m_resize(Y,s,A->n);
  4198. -   MEM_STAT_REG(Y,TYPE_MAT);
  4199. -   /* pointers to rows of Y */
  4200. -   y0.dim = y0.max_dim = A->n;
  4201. -   y1.dim = y1.max_dim = A->n;
  4202. -
  4203. -   m_zero(Y);
  4204. -   m_zero(out);
  4205. -   
  4206. -#define Z(k)     ((k) & 1 ? tmp : &y0)
  4207. -#define ZZ(k)    ((k) & 1 ? tmp->ve : y0.ve)
  4208. -
  4209. -   for( j = 0; j < A->n; j++)
  4210. -   {
  4211. -      if( j > 0 )
  4212. -    Y->me[0][j-1] = 0.0;
  4213. -      Y->me[0][j] = 1.0;
  4214. -
  4215. -      y0.ve = Y->me[0];
  4216. -      for (k = 0; k < s-1; k++)
  4217. -      {
  4218. -     y1.ve = Y->me[k+1];
  4219. -     mv_mlt(A,&y0,&y1);
  4220. -     y0.ve = y1.ve;
  4221. -      }
  4222. -      
  4223. -      y0.ve = out->me[j];
  4224. -
  4225. -      t = s*r;
  4226. -      for ( l = 0; l <= q-t; l++ )
  4227. -    __mltadd__(y0.ve,Y->me[l],a->ve[t+l],Y->n);
  4228. -      
  4229. -      for (k=1; k <= r; k++)
  4230. -      {
  4231. -     mv_mlt(Apow,Z(k-1),Z(k)); 
  4232. -     t = s*(r-k);
  4233. -     for (l=0; l < s; l++)
  4234. -       __mltadd__(ZZ(k),Y->me[l],a->ve[t+l],Y->n);
  4235. -      }
  4236. -      if (Z(k) == &y0) v_copy(tmp,&y0);
  4237. -   }
  4238. -
  4239. -   m_transp(out,out);
  4240. -   
  4241. -   return out;
  4242. -}
  4243. -
  4244. -
  4245. //GO.SYSIN DD mfunc.c
  4246. echo bdfactor.c 1>&2
  4247. sed >bdfactor.c <<'//GO.SYSIN DD bdfactor.c' 's/^-//'
  4248. -
  4249. -/**************************************************************************
  4250. -**
  4251. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4252. -**
  4253. -**                 Meschach Library
  4254. -** 
  4255. -** This Meschach Library is provided "as is" without any express 
  4256. -** or implied warranty of any kind with respect to this software. 
  4257. -** In particular the authors shall not be liable for any direct, 
  4258. -** indirect, special, incidental or consequential damages arising 
  4259. -** in any way from use of the software.
  4260. -** 
  4261. -** Everyone is granted permission to copy, modify and redistribute this
  4262. -** Meschach Library, provided:
  4263. -**  1.  All copies contain this copyright notice.
  4264. -**  2.  All modified copies shall carry a notice stating who
  4265. -**      made the last modification and the date of such modification.
  4266. -**  3.  No charge is made for this software or works derived from it.  
  4267. -**      This clause shall not be construed as constraining other software
  4268. -**      distributed on the same medium as this software, nor is a
  4269. -**      distribution fee considered a charge.
  4270. -**
  4271. -***************************************************************************/
  4272. -
  4273. -
  4274. -/*
  4275. -  Band matrix factorisation routines
  4276. -  */
  4277. -
  4278. -/* bdfactor.c  18/11/93 */
  4279. -static    char    rcsid[] = "$Id: bdfactor.c,v 1.5 1994/05/17 23:18:32 des Exp $";
  4280. -
  4281. -#include    <stdio.h>
  4282. -#include    <math.h>
  4283. -#include        "matrix2.h"
  4284. -
  4285. -
  4286. -/* generate band matrix 
  4287. -   for a matrix  with n columns,
  4288. -   lb subdiagonals and ub superdiagonals;
  4289. -
  4290. -   Way of saving a band of a matrix:
  4291. -   first we save subdiagonals (from 0 to lb-1);
  4292. -   then main diagonal (in the lb row)
  4293. -   and then superdiagonals (from lb+1 to lb+ub)
  4294. -   in such a way that the elements which were previously
  4295. -   in one column are now also in one column
  4296. -*/
  4297. -
  4298. -BAND *bd_get(lb,ub,n)
  4299. -int lb, ub, n;
  4300. -{
  4301. -   BAND *A;
  4302. -
  4303. -   if (lb < 0 || ub < 0 || n <= 0)
  4304. -     error(E_NEG,"bd_get");
  4305. -
  4306. -   if ((A = NEW(BAND)) == (BAND *)NULL)
  4307. -     error(E_MEM,"bd_get");
  4308. -   else if (mem_info_is_on()) {
  4309. -      mem_bytes(TYPE_BAND,0,sizeof(BAND));
  4310. -      mem_numvar(TYPE_BAND,1);
  4311. -   }
  4312. -
  4313. -   lb = A->lb = min(n-1,lb);
  4314. -   ub = A->ub = min(n-1,ub);
  4315. -   A->mat = m_get(lb+ub+1,n);
  4316. -   return A;
  4317. -}
  4318. -
  4319. -int bd_free(A)
  4320. -BAND *A;
  4321. -{
  4322. -   if ( A == (BAND *)NULL || A->lb < 0 || A->ub < 0 )
  4323. -     /* don't trust it */
  4324. -     return (-1);
  4325. -
  4326. -   if (A->mat) m_free(A->mat);
  4327. -
  4328. -   if (mem_info_is_on()) {
  4329. -      mem_bytes(TYPE_BAND,sizeof(BAND),0);
  4330. -      mem_numvar(TYPE_BAND,-1);
  4331. -   }
  4332. -
  4333. -   free((char *)A);
  4334. -   return 0;
  4335. -}
  4336. -
  4337. -
  4338. -/* resize band matrix */
  4339. -
  4340. -BAND *bd_resize(A,new_lb,new_ub,new_n)
  4341. -BAND *A;
  4342. -int new_lb,new_ub,new_n;
  4343. -{
  4344. -   int lb,ub,i,j,l,shift,umin;
  4345. -   Real **Av;
  4346. -
  4347. -   if (new_lb < 0 || new_ub < 0 || new_n <= 0)
  4348. -     error(E_NEG,"bd_resize");
  4349. -   if ( ! A )
  4350. -     return bd_get(new_lb,new_ub,new_n);
  4351. -    if ( A->lb+A->ub+1 > A->mat->m )
  4352. -    error(E_INTERN,"bd_resize");
  4353. -
  4354. -   if ( A->lb == new_lb && A->ub == new_ub && A->mat->n == new_n )
  4355. -    return A;
  4356. -
  4357. -   lb = A->lb;
  4358. -   ub = A->ub;
  4359. -   Av = A->mat->me;
  4360. -   umin = min(ub,new_ub);
  4361. -
  4362. -    /* ensure that unused triangles at edges are zero'd */
  4363. -
  4364. -   for ( i = 0; i < lb; i++ )
  4365. -      for ( j = A->mat->n - lb + i; j < A->mat->n; j++ )
  4366. -    Av[i][j] = 0.0;  
  4367. -    for ( i = lb+1,l=1; l <= umin; i++,l++ )
  4368. -      for ( j = 0; j < l; j++ )
  4369. -    Av[i][j] = 0.0; 
  4370. -
  4371. -   new_lb = A->lb = min(new_lb,new_n-1);
  4372. -   new_ub = A->ub = min(new_ub,new_n-1);
  4373. -   A->mat = m_resize(A->mat,new_lb+new_ub+1,new_n);
  4374. -   Av = A->mat->me;
  4375. -
  4376. -   /* if new_lb != lb then move the rows to get the main diag 
  4377. -      in the new_lb row */
  4378. -
  4379. -   if (new_lb > lb) {
  4380. -      shift = new_lb-lb;
  4381. -
  4382. -      for (i=lb+umin, l=i+shift; i >= 0; i--,l--)
  4383. -    MEM_COPY(Av[i],Av[l],new_n*sizeof(Real));
  4384. -      for (l=shift-1; l >= 0; l--)
  4385. -    __zero__(Av[l],new_n);
  4386. -   }
  4387. -   else { /* new_lb < lb */
  4388. -      shift = lb - new_lb;
  4389. -
  4390. -      for (i=shift, l=0; i <= lb+umin; i++,l++)
  4391. -    MEM_COPY(Av[i],Av[l],new_n*sizeof(Real));
  4392. -      for (i=lb+umin+1; i <= new_lb+new_ub; i++)
  4393. -    __zero__(Av[i],new_n);
  4394. -   }
  4395. -
  4396. -   return A;
  4397. -}
  4398. -
  4399. -
  4400. -
  4401. -BAND *bd_copy(A,B)
  4402. -BAND *A,*B;
  4403. -{
  4404. -   int lb,ub,i,j,n;
  4405. -   
  4406. -   if ( !A )
  4407. -     error(E_NULL,"bd_copy");
  4408. -
  4409. -   if (A == B) return B;
  4410. -   
  4411. -   n = A->mat->n;
  4412. -   if ( !B )
  4413. -     B = bd_get(A->lb,A->ub,n);
  4414. -   else if (B->lb != A->lb || B->ub != A->ub || B->mat->n != n )
  4415. -     B = bd_resize(B,A->lb,A->ub,n);
  4416. -   
  4417. -   if (A->mat == B->mat) return B;
  4418. -   ub = B->ub = A->ub;
  4419. -   lb = B->lb = A->lb;
  4420. -
  4421. -   for ( i=0, j=n-lb; i <= lb; i++, j++ )
  4422. -     MEM_COPY(A->mat->me[i],B->mat->me[i],j*sizeof(Real));   
  4423. -
  4424. -   for ( i=lb+1, j=1; i <= lb+ub; i++, j++ )
  4425. -     MEM_COPY(A->mat->me[i]+j,B->mat->me[i]+j,(n - j)*sizeof(Real));     
  4426. -
  4427. -   return B;
  4428. -}
  4429. -
  4430. -
  4431. -/* copy band matrix to a square matrix */
  4432. -MAT *band2mat(bA,A)
  4433. -BAND *bA;
  4434. -MAT *A;
  4435. -{
  4436. -   int i,j,l,n,n1;
  4437. -   int lb, ub;
  4438. -   Real **bmat;
  4439. -
  4440. -   if ( !bA || !A)
  4441. -     error(E_NULL,"band2mat");
  4442. -   if ( bA->mat == A )
  4443. -     error(E_INSITU,"band2mat");
  4444. -
  4445. -   ub = bA->ub;
  4446. -   lb = bA->lb;
  4447. -   n = bA->mat->n;
  4448. -   n1 = n-1;
  4449. -   bmat = bA->mat->me;
  4450. -
  4451. -   A = m_resize(A,n,n);
  4452. -   m_zero(A);
  4453. -
  4454. -   for (j=0; j < n; j++)
  4455. -     for (i=min(n1,j+lb),l=lb+j-i; i >= max(0,j-ub); i--,l++)
  4456. -       A->me[i][j] = bmat[l][j];
  4457. -
  4458. -   return A;
  4459. -}
  4460. -
  4461. -/* copy a square matrix to a band matrix with 
  4462. -   lb subdiagonals and ub superdiagonals */
  4463. -BAND *mat2band(A,lb,ub,bA)
  4464. -BAND *bA;
  4465. -MAT *A;
  4466. -int lb, ub;
  4467. -{
  4468. -   int i, j, l, n1;
  4469. -   Real **bmat;
  4470. -   
  4471. -   if (! A || ! bA)
  4472. -     error(E_NULL,"mat2band");
  4473. -   if (ub < 0 || lb < 0)
  4474. -     error(E_SIZES,"mat2band");
  4475. -   if (bA->mat == A)
  4476. -     error(E_INSITU,"mat2band");
  4477. -
  4478. -   n1 = A->n-1;
  4479. -   lb = min(n1,lb);
  4480. -   ub = min(n1,ub);
  4481. -   bA = bd_resize(bA,lb,ub,n1+1);
  4482. -   bmat = bA->mat->me;
  4483. -
  4484. -   for (j=0; j <= n1; j++)
  4485. -     for (i=min(n1,j+lb),l=lb+j-i; i >= max(0,j-ub); i--,l++)
  4486. -       bmat[l][j] = A->me[i][j];
  4487. -
  4488. -   return bA;
  4489. -}
  4490. -
  4491. -
  4492. -
  4493. -/* transposition of matrix in;
  4494. -   out - matrix after transposition;
  4495. -   can be done in situ
  4496. -*/
  4497. -
  4498. -BAND *bd_transp(in,out)
  4499. -BAND *in, *out;
  4500. -{
  4501. -   int i, j, jj, l, k, lb, ub, lub, n, n1;
  4502. -   int in_situ;
  4503. -   Real  **in_v, **out_v;
  4504. -   
  4505. -   if ( in == (BAND *)NULL || in->mat == (MAT *)NULL )
  4506. -     error(E_NULL,"bd_transp");
  4507. -
  4508. -   lb = in->lb;
  4509. -   ub = in->ub;
  4510. -   lub = lb+ub;
  4511. -   n = in->mat->n;
  4512. -   n1 = n-1;
  4513. -
  4514. -   in_situ = ( in == out );
  4515. -   if ( ! in_situ )
  4516. -       out = bd_resize(out,ub,lb,n);
  4517. -   else
  4518. -   {   /* only need to swap lb and ub fields */
  4519. -       out->lb = ub;
  4520. -       out->ub = lb;
  4521. -   }
  4522. -
  4523. -   in_v = in->mat->me;
  4524. -   
  4525. -   if (! in_situ) {
  4526. -      int sh_in,sh_out; 
  4527. -
  4528. -      out_v = out->mat->me;
  4529. -      for (i=0, l=lub, k=lb-i; i <= lub; i++,l--,k--) {
  4530. -     sh_in = max(-k,0);
  4531. -     sh_out = max(k,0);
  4532. -     MEM_COPY(&(in_v[i][sh_in]),&(out_v[l][sh_out]),
  4533. -          (n-sh_in-sh_out)*sizeof(Real));
  4534. -     /**********************************
  4535. -     for (j=n1-sh_out, jj=n1-sh_in; j >= sh_in; j--,jj--) {
  4536. -        out_v[l][jj] = in_v[i][j];
  4537. -     }
  4538. -     **********************************/
  4539. -      }
  4540. -   }
  4541. -   else if (ub == lb) {
  4542. -      Real tmp;
  4543. -
  4544. -      for (i=0, l=lub, k=lb-i; i < lb; i++,l--,k--) {
  4545. -     for (j=n1-k, jj=n1; j >= 0; j--,jj--) {
  4546. -        tmp = in_v[l][jj];
  4547. -        in_v[l][jj] = in_v[i][j];
  4548. -        in_v[i][j] = tmp;
  4549. -     }
  4550. -      }
  4551. -   }
  4552. -   else if (ub > lb) {  /* hence i-ub <= 0 & l-lb >= 0 */
  4553. -      int p,pp,lbi;
  4554. -      
  4555. -      for (i=0, l=lub; i < (lub+1)/2; i++,l--) {
  4556. -     lbi = lb-i;
  4557. -     for (j=l-lb, jj=0, p=max(-lbi,0), pp = max(l-ub,0); j <= n1; 
  4558. -          j++,jj++,p++,pp++) {
  4559. -        in_v[l][pp] = in_v[i][p];
  4560. -        in_v[i][jj] = in_v[l][j];
  4561. -     }
  4562. -     for (  ; p <= n1-max(lbi,0); p++,pp++)
  4563. -       in_v[l][pp] = in_v[i][p];
  4564. -      }
  4565. -      
  4566. -      if (lub%2 == 0) { /* shift only */
  4567. -     i = lub/2;
  4568. -     for (j=max(i-lb,0), jj=0; jj <= n1-ub+i; j++,jj++) 
  4569. -       in_v[i][jj] = in_v[i][j];
  4570. -      }
  4571. -   }
  4572. -   else {      /* ub < lb, hence ub-l <= 0 & lb-i >= 0 */
  4573. -      int p,pp,ubi;
  4574. -
  4575. -      for (i=0, l=lub; i < (lub+1)/2; i++,l--) {
  4576. -     ubi = i-ub;
  4577. -     for (j=n1-max(lb-l,0), jj=n1-max(-ubi,0), p=n1-lb+i, pp=n1;
  4578. -          p >= 0; j--, jj--, pp--, p--) {
  4579. -        in_v[i][jj] = in_v[l][j];
  4580. -        in_v[l][pp] = in_v[i][p];
  4581. -     }
  4582. -     for (  ; jj >= max(ubi,0); j--, jj--)
  4583. -       in_v[i][jj] = in_v[l][j];
  4584. -      }
  4585. -
  4586. -      if (lub%2 == 0) {  /* shift only */
  4587. -     i = lub/2;
  4588. -     for (j=n1-lb+i, jj=n1-max(ub-i,0); j >= 0; j--, jj--) 
  4589. -        in_v[i][jj] = in_v[i][j];
  4590. -      }
  4591. -   }
  4592. -
  4593. -   return out;
  4594. -}
  4595. -
  4596. -
  4597. -
  4598. -/* bdLUfactor -- gaussian elimination with partial pivoting
  4599. -   -- on entry, the matrix A in band storage with elements 
  4600. -      in rows 0 to lb+ub; 
  4601. -      The jth column of A is stored in the jth column of 
  4602. -      band A (bA) as follows:
  4603. -      bA->mat->me[lb+j-i][j] = A->me[i][j] for 
  4604. -      max(0,j-lb) <= i <= min(A->n-1,j+ub);
  4605. -   -- on exit: U is stored as an upper triangular matrix
  4606. -      with lb+ub superdiagonals in rows lb to 2*lb+ub, 
  4607. -      and the matrix L is stored in rows 0 to lb-1.
  4608. -      Matrix U is permuted, whereas L is not permuted !!!
  4609. -      Therefore we save some memory.
  4610. -   */
  4611. -BAND    *bdLUfactor(bA,pivot)
  4612. -BAND    *bA;
  4613. -PERM    *pivot;
  4614. -{
  4615. -   int    i, j, k, l, n, n1, lb, ub, lub, k_end, k_lub;
  4616. -   int    i_max, shift;
  4617. -   Real    **bA_v;
  4618. -   Real max1, temp;
  4619. -   
  4620. -   if ( bA==(BAND *)NULL || pivot==(PERM *)NULL )
  4621. -     error(E_NULL,"bdLUfactor");
  4622. -
  4623. -   lb = bA->lb;
  4624. -   ub = bA->ub;
  4625. -   lub = lb+ub;
  4626. -   n = bA->mat->n;
  4627. -   n1 = n-1;
  4628. -   lub = lb+ub;
  4629. -
  4630. -   if ( pivot->size != n )
  4631. -     error(E_SIZES,"bdLUfactor");
  4632. -
  4633. -   
  4634. -   /* initialise pivot with identity permutation */
  4635. -   for ( i=0; i < n; i++ )
  4636. -     pivot->pe[i] = i;
  4637. -
  4638. -   /* extend band matrix */
  4639. -   /* extended part is filled with zeros */
  4640. -   bA = bd_resize(bA,lb,min(n1,lub),n);
  4641. -   bA_v = bA->mat->me;
  4642. -
  4643. -
  4644. -   /* main loop */
  4645. -
  4646. -   for ( k=0; k < n1; k++ )
  4647. -   {
  4648. -      k_end = max(0,lb+k-n1);
  4649. -      k_lub = min(k+lub,n1);
  4650. -
  4651. -      /* find the best pivot row */
  4652. -      
  4653. -      max1 = 0.0;    
  4654. -      i_max = -1;
  4655. -      for ( i=lb; i >= k_end; i-- ) {
  4656. -     temp = fabs(bA_v[i][k]);
  4657. -     if ( temp > max1 )
  4658. -     { max1 = temp;    i_max = i; }
  4659. -      }
  4660. -      
  4661. -      /* if no pivot then ignore column k... */
  4662. -      if ( i_max == -1 )
  4663. -    continue;
  4664. -      
  4665. -      /* do we pivot ? */
  4666. -      if ( i_max != lb )    /* yes we do... */
  4667. -      {
  4668. -     /* save transposition using non-shifted indices */
  4669. -     shift = lb-i_max;
  4670. -     px_transp(pivot,k+shift,k);
  4671. -     for ( i=lb, j=k; j <= k_lub; i++,j++ )
  4672. -     {
  4673. -        temp = bA_v[i][j];
  4674. -        bA_v[i][j] = bA_v[i-shift][j];
  4675. -        bA_v[i-shift][j] = temp;
  4676. -     }
  4677. -      }
  4678. -      
  4679. -      /* row operations */
  4680. -      for ( i=lb-1; i >= k_end; i-- ) {
  4681. -     temp = bA_v[i][k] /= bA_v[lb][k];
  4682. -     shift = lb-i;
  4683. -     for ( j=k+1,l=i+1; j <= k_lub; l++,j++ )
  4684. -       bA_v[l][j] -= temp*bA_v[l+shift][j];
  4685. -      }
  4686. -   }
  4687. -   
  4688. -   return bA;
  4689. -}
  4690. -
  4691. -
  4692. -/* bdLUsolve -- given an LU factorisation in bA, solve bA*x=b */
  4693. -/* pivot is changed upon return  */
  4694. -VEC    *bdLUsolve(bA,pivot,b,x)
  4695. -BAND    *bA;
  4696. -PERM    *pivot;
  4697. -VEC    *b,*x;
  4698. -{
  4699. -   int i,j,l,n,n1,pi,lb,ub,jmin, maxj;
  4700. -   Real c;
  4701. -   Real **bA_v;
  4702. -
  4703. -   if ( bA==(BAND *)NULL || b==(VEC *)NULL || pivot==(PERM *)NULL )
  4704. -     error(E_NULL,"bdLUsolve");
  4705. -   if ( bA->mat->n != b->dim || bA->mat->n != pivot->size)
  4706. -     error(E_SIZES,"bdLUsolve");
  4707. - 
  4708. -   lb = bA->lb;
  4709. -   ub = bA->ub;
  4710. -   n = b->dim;
  4711. -   n1 = n-1;
  4712. -   bA_v = bA->mat->me;
  4713. -
  4714. -   x = v_resize(x,b->dim);
  4715. -   px_vec(pivot,b,x);
  4716. -
  4717. -   /* solve Lx = b; implicit diagonal = 1 
  4718. -      L is not permuted, therefore it must be permuted now
  4719. -    */
  4720. -   
  4721. -   px_inv(pivot,pivot);
  4722. -   for (j=0; j < n; j++) {
  4723. -      jmin = j+1;
  4724. -      c = x->ve[j];
  4725. -      maxj = max(0,j+lb-n1);
  4726. -      for (i=jmin,l=lb-1; l >= maxj; i++,l--) {
  4727. -     if ( (pi = pivot->pe[i]) < jmin) 
  4728. -       pi = pivot->pe[i] = pivot->pe[pi];
  4729. -     x->ve[pi] -= bA_v[l][j]*c;
  4730. -      }
  4731. -   }
  4732. -
  4733. -   /* solve Ux = b; explicit diagonal */
  4734. -
  4735. -   x->ve[n1] /= bA_v[lb][n1];
  4736. -   for (i=n-2; i >= 0; i--) {
  4737. -      c = x->ve[i];
  4738. -      for (j=min(n1,i+ub), l=lb+j-i; j > i; j--,l--)
  4739. -    c -= bA_v[l][j]*x->ve[j];
  4740. -      x->ve[i] = c/bA_v[lb][i];
  4741. -   }
  4742. -   
  4743. -   return (x);
  4744. -}
  4745. -
  4746. -/* LDLfactor -- L.D.L' factorisation of A in-situ;
  4747. -   A is a band matrix
  4748. -   it works using only lower bandwidth & main diagonal
  4749. -   so it is possible to set A->ub = 0
  4750. - */
  4751. -
  4752. -BAND *bdLDLfactor(A)
  4753. -BAND *A;
  4754. -{
  4755. -   int i,j,k,n,n1,lb,ki,jk,ji,lbkm,lbkp;
  4756. -   Real **Av;
  4757. -   Real c, cc;
  4758. -
  4759. -   if ( ! A )
  4760. -     error(E_NULL,"bdLDLfactor");
  4761. -
  4762. -   if (A->lb == 0) return A;
  4763. -
  4764. -   lb = A->lb;
  4765. -   n = A->mat->n;
  4766. -   n1 = n-1;
  4767. -   Av = A->mat->me;
  4768. -   
  4769. -   for (k=0; k < n; k++) {    
  4770. -      lbkm = lb-k;
  4771. -      lbkp = lb+k;
  4772. -
  4773. -      /* matrix D */
  4774. -      c = Av[lb][k];
  4775. -      for (j=max(0,-lbkm), jk=lbkm+j; j < k; j++, jk++) {
  4776. -     cc = Av[jk][j];
  4777. -     c -= Av[lb][j]*cc*cc;
  4778. -      }
  4779. -      if (c == 0.0)
  4780. -    error(E_SING,"bdLDLfactor");
  4781. -      Av[lb][k] = c;
  4782. -
  4783. -      /* matrix L */
  4784. -      
  4785. -      for (i=min(n1,lbkp), ki=lbkp-i; i > k; i--,ki++) {
  4786. -     c = Av[ki][k];
  4787. -     for (j=max(0,i-lb), ji=lb+j-i, jk=lbkm+j; j < k;
  4788. -          j++, ji++, jk++)
  4789. -       c -= Av[lb][j]*Av[ji][j]*Av[jk][j];
  4790. -     Av[ki][k] = c/Av[lb][k];
  4791. -      }
  4792. -   }
  4793. -   
  4794. -   return A;
  4795. -}
  4796. -
  4797. -/* solve A*x = b, where A is factorized by 
  4798. -   Choleski LDL^T factorization */
  4799. -VEC    *bdLDLsolve(A,b,x)
  4800. -BAND   *A;
  4801. -VEC    *b, *x;
  4802. -{
  4803. -   int i,j,l,n,n1,lb,ilb;
  4804. -   Real **Av, *Avlb;
  4805. -   Real c;
  4806. -
  4807. -   if ( ! A || ! b )
  4808. -     error(E_NULL,"bdLDLsolve");
  4809. -   if ( A->mat->n != b->dim )
  4810. -     error(E_SIZES,"bdLDLsolve");
  4811. -
  4812. -   n = A->mat->n;
  4813. -   n1 = n-1;
  4814. -   x = v_resize(x,n);
  4815. -   lb = A->lb;
  4816. -   Av = A->mat->me;  
  4817. -   Avlb = Av[lb];
  4818. -   
  4819. -   /* solve L*y = b */
  4820. -   x->ve[0] = b->ve[0];
  4821. -   for (i=1; i < n; i++) {
  4822. -      ilb = i-lb;
  4823. -      c = b->ve[i];
  4824. -      for (j=max(0,ilb), l=j-ilb; j < i; j++,l++)
  4825. -    c -= Av[l][j]*x->ve[j];
  4826. -      x->ve[i] = c;
  4827. -   }
  4828. -
  4829. -   /* solve D*z = y */
  4830. -   for (i=0; i < n; i++) 
  4831. -     x->ve[i] /= Avlb[i];
  4832. -
  4833. -   /* solve L^T*x = z */
  4834. -   for (i=n-2; i >= 0; i--) {
  4835. -      ilb = i+lb;
  4836. -      c = x->ve[i];
  4837. -      for (j=min(n1,ilb), l=ilb-j; j > i; j--,l++)
  4838. -    c -= Av[l][i]*x->ve[j];
  4839. -      x->ve[i] = c;
  4840. -   }
  4841. -
  4842. -   return x;
  4843. -}
  4844. -
  4845. -
  4846. -
  4847. //GO.SYSIN DD bdfactor.c
  4848.